{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "85a8fc3f",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# 导言区"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f397dd62",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.constants as sc\n",
    "from CoolProp.CoolProp import PropsSI as psi\n",
    "import sympy as sp\n",
    "from sympy import symbols, diff, Function, dsolve, solve, pi, integrate\n",
    "from scipy.integrate import quad\n",
    "import matplotlib.pyplot as plt\n",
    "from Appendix import Appendix5_air_physical_properties as ap5\n",
    "from Appendix.Appendix4_lambda_ import get_lambda_\n",
    "from Functions.SteadyStateConduction import *\n",
    "from Functions.UnsteadyStateConduction import *\n",
    "from Appendix import Appendix4_lambda_ as ap4\n",
    "from scipy.integrate import solve_bvp\n",
    "from scipy.special import iv\n",
    "from Functions.Self_defined import find_nearest"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79c7c6af",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 习题02-01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "daa3b372",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t = 238.20 C\n"
     ]
    }
   ],
   "source": [
    "t_boiler = 111\n",
    "q = 42400\n",
    "delta = 3e-3\n",
    "lambda_ = 1\n",
    "\n",
    "t = root(lambda t: q - lambda_ * (t - t_boiler) / delta,\n",
    "         t_boiler + 10).x[0]\n",
    "print(f't = {t:.2f} C')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ea4140c",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 习题02-02"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "89fb1fcc",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Q = 1285717.03 J\n"
     ]
    }
   ],
   "source": [
    "delta_1, delta_2, delta_3 = 0.794e-3, 152e-3, 9.5e-3\n",
    "lambda_1, lambda_2, lambda_3 = 45, 0.07, 0.1\n",
    "area = 37.2\n",
    "t_in = -2\n",
    "t_out = 30\n",
    "h_in = 1.5\n",
    "h_out = 2.5\n",
    "\n",
    "tau = 1*3600\n",
    "\n",
    "R_1 = delta_1 / (lambda_1 * area)\n",
    "R_2 = delta_2 / (lambda_2 * area)\n",
    "R_3 = delta_3 / (lambda_3 * area)\n",
    "R_in = 1 / (h_in * area)\n",
    "R_out = 1 / (h_out * area)\n",
    "R_total = R_1 + R_2 + R_3 + R_in + R_out\n",
    "\n",
    "q = (t_out - t_in) / R_total\n",
    "Q = q * tau\n",
    "print(f'Q = {Q:.2f} J')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65d7b32d",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 习题02-03"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4094abc2",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "保温层的厚度为 36 mm\n"
     ]
    }
   ],
   "source": [
    "delta = 20e-3\n",
    "lambda_ = 1.3\n",
    "q = 1500\n",
    "lambda_insu = 0.08\n",
    "t_in = 750\n",
    "t_out = 55\n",
    "\n",
    "R_1 = delta / lambda_\n",
    "R_total = (t_in - t_out) / q\n",
    "R_2 = R_total - R_1\n",
    "delta_insu = lambda_insu * R_2\n",
    "print(f'保温层的厚度为 {delta_insu*1000:.0f} mm')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1630f29",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 习题02-04"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "cbee7d1c",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "delta_A = 0.065 m\n",
      "delta_B = 0.032 m\n"
     ]
    }
   ],
   "source": [
    "delta_ratio_A_B = 2\n",
    "lambda_A = 0.08\n",
    "lambda_B = 0.05\n",
    "t_f1 = 400\n",
    "h_1 = 50\n",
    "t_w = 50\n",
    "t_f2 = 25\n",
    "h_2 = 9.5\n",
    "\n",
    "q = h_2 * (t_w - t_f2)\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    delta_A, delta_B, R_total = p\n",
    "    expr1 = delta_A - delta_B * delta_ratio_A_B\n",
    "    expr2 = R_total - 1/h_1 - 1/h_2 - delta_A/lambda_A - delta_B/lambda_B\n",
    "    expr3 = q - (t_f1 - t_f2) / R_total\n",
    "    return expr1, expr2, expr3\n",
    "\n",
    "\n",
    "guess_value = (0.02, 0.01, 1)\n",
    "delta_A, delta_B, R_total = root(expressions, guess_value).x\n",
    "print(f'delta_A = {delta_A:.3f} m')\n",
    "print(f'delta_B = {delta_B:.3f} m')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61e0b275",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 习题02-06"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "deb56db1",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "该燃烧室工作于安全温度范围内\n"
     ]
    }
   ],
   "source": [
    "d_o = 130e-3\n",
    "delta = 2.1e-3\n",
    "lambda_ = 23.2\n",
    "t_o = 240\n",
    "q = 4.8e6\n",
    "t_i_max = 700\n",
    "\n",
    "delta_t = q * delta / lambda_\n",
    "t_i = t_o + delta_t\n",
    "if t_i < t_i_max:\n",
    "    print('该燃烧室工作于安全温度范围内')\n",
    "else:\n",
    "    print('该燃烧室工作于安全温度范围外')"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-07"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t_bot = 105.96 C\n",
      "t_top = 105.82 C\n"
     ]
    }
   ],
   "source": [
    "P_total = 1000\n",
    "eta = 85/100\n",
    "delta = 3e-3\n",
    "d = 200e-3\n",
    "lambda_ = 18\n",
    "h = 2500\n",
    "t_f = 95\n",
    "\n",
    "P = P_total * eta\n",
    "area = np.pi * d**2 / 4\n",
    "R_1 = 1 / (h * area)\n",
    "R_2 = delta / lambda_\n",
    "R = R_1 + R_2\n",
    "t_bot = t_f + P * R\n",
    "t_top = t_f + P * R_1\n",
    "print(f't_bot = {t_bot:.2f} C')\n",
    "print(f't_top = {t_top:.2f} C')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-08"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lambda_ = 1.43770613178973*lambda_1\n",
      "标准材料的导热系数随温度发生变化时，Delta_t1与Delta_t2不等\n"
     ]
    }
   ],
   "source": [
    "t_h = 400\n",
    "t_c = 300\n",
    "Delta_t = 2.49\n",
    "Delta_t1 = 3.56\n",
    "Delta_t2 = 3.60\n",
    "\n",
    "q, A, l, lambda_, lambda_1, lambda_2 = symbols('q A l lambda_ lambda_1 lambda_2')\n",
    "\n",
    "xpr1 = q - Delta_t1 * A * lambda_1 / l\n",
    "xpr2 = q - Delta_t2 * A * lambda_2 / l\n",
    "xpr3 = q - Delta_t * A * lambda_ / l\n",
    "\n",
    "result = solve([xpr1, xpr2, xpr3], [lambda_1, lambda_, lambda_2])\n",
    "lambda_ = 2 * result[lambda_] / (result[lambda_1] + result[lambda_2]) * lambda_1\n",
    "print(f'lambda_ = {lambda_}')\n",
    "\n",
    "if Delta_t1 == Delta_t2:\n",
    "    print('标准材料的导热系数随温度发生变化时，Delta_t1与Delta_t2相等')\n",
    "else:\n",
    "    print('标准材料的导热系数随温度发生变化时，Delta_t1与Delta_t2不等')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-09"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "热损失为41.95 W\n",
      "其热损失是双层玻璃的44.62倍\n"
     ]
    }
   ],
   "source": [
    "delta_1 = delta_2 = 6e-3\n",
    "delta_air = 8e-3\n",
    "t_in = 20\n",
    "t_out = -20\n",
    "width = length = 60e-2\n",
    "lambda_ = 0.78\n",
    "\n",
    "t = (t_in + t_out) / 2\n",
    "lambda_air = np.interp(t, ap5.temperature_list, ap5.lambda_list)\n",
    "# 也可以利用CoolProp计算lambda的值\n",
    "# lambda_air2 = psi('L', 'T', t+273.15, 'P', 101325, 'Air')\n",
    "\n",
    "area = width * length\n",
    "R_1 = delta_1 / (lambda_ * area)\n",
    "R_2 = delta_2 / (lambda_ * area)\n",
    "R_air = delta_air / (lambda_air * area)\n",
    "\n",
    "R = R_1 + R_2 + R_air\n",
    "q = (t_in - t_out) / R\n",
    "print(f'热损失为{q:.2f} W')\n",
    "\n",
    "R_p = delta_1 / (lambda_ * area)\n",
    "q_p = (t_in - t_out) / R_p\n",
    "times = q_p / q\n",
    "print(f'其热损失是双层玻璃的{times:.2f}倍')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-10"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "q = 49.99 W/m^2\n"
     ]
    }
   ],
   "source": [
    "delta_g = 3e-3\n",
    "delta_air = 6e-3\n",
    "lambda_g = 0.8\n",
    "t_i = 15\n",
    "t_o = -10\n",
    "\n",
    "air_pressure = sc.atm\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    t_12, t_23, t_34, t_45, q, lambda_air1, lambda_air2 = p\n",
    "    xpr1 = q - (t_i - t_12) * lambda_g / delta_g\n",
    "    xpr2 = q - (t_12 - t_23) * lambda_air1 / delta_air\n",
    "    xpr3 = q - (t_23 - t_34) * lambda_g / delta_g\n",
    "    xpr4 = q - (t_34 - t_45) * lambda_air2 / delta_air\n",
    "    xpr5 = q - (t_45 - t_o) * lambda_g / delta_g\n",
    "    xpr6 = lambda_air1 - psi('L', 'T', sc.convert_temperature(((t_12 + t_23) / 2), 'C', 'K'), 'P', air_pressure, 'Air')\n",
    "    xpr7 = lambda_air2 - psi('L', 'T', sc.convert_temperature(((t_34 + t_45) / 2), 'C', 'K'), 'P', air_pressure, 'Air')\n",
    "    return [xpr1, xpr2, xpr3, xpr4, xpr5, xpr6, xpr7]\n",
    "\n",
    "\n",
    "guess_value = [t_i- 1, t_i-2, t_i-3, t_i-4, 1, lambda_g, lambda_g]\n",
    "t_12, t_23, t_34, t_45, q, lambda_air1, lambda_air2 = root(expressions, guess_value).x\n",
    "print(f'q = {q:.2f} W/m^2')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "id": "35728adb",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 习题02-11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "合金的最高温度为 1078.39 C。\n",
      "合金可以安全地工作\n"
     ]
    }
   ],
   "source": [
    "lambda_ceramics = 1.3\n",
    "t_max = 1250\n",
    "lambda_alloy = 25\n",
    "R_c = 1e-4  # m^2-K/W\n",
    "t_gas = 1700\n",
    "h_gas = 1000\n",
    "t_air = 400\n",
    "h_air = 500\n",
    "delta_ceramics = 1e-3   # 缺失\n",
    "delta_alloy = 1e-3  # 缺失\n",
    "\n",
    "R_1 = 1 / h_gas\n",
    "R_2 = delta_ceramics / lambda_ceramics  # 缺少条件：陶瓷层的厚度\n",
    "R_3 = delta_alloy / lambda_alloy    # 缺少条件：合金层的厚度\n",
    "R_4 = 1 / h_air\n",
    "R_total = R_1 + R_2 + R_c + R_3 + R_4\n",
    "\n",
    "q = (t_gas - t_air) / R_total\n",
    "t_alloy_max = t_air + q * (R_3 + R_4)\n",
    "print(f'合金的最高温度为 {t_alloy_max:.2f} C。')\n",
    "\n",
    "if t_alloy_max > t_max:\n",
    "    print('合金不能安全地工作')\n",
    "else:\n",
    "    print('合金可以安全地工作')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-12"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "q = 2942.86 W/m^2\n"
     ]
    }
   ],
   "source": [
    "delta_s = 1e-3  # 基板的厚度\n",
    "delta_f = 0.2e-3    # 薄膜的厚度\n",
    "t_f = 20\n",
    "h = 40\n",
    "t_1 = 30\n",
    "t_0 = 60\n",
    "lambda_f = 0.02\n",
    "lambda_s = 0.06\n",
    "\n",
    "R_1 = 1 / h\n",
    "R_2 = delta_f / lambda_f\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    q_1, t_2 = p\n",
    "    xpr1 = q_1 - h * (t_2 - t_f)\n",
    "    xpr2 = q_1 - (t_0 - t_2) * lambda_f / delta_f\n",
    "    return [xpr1, xpr2]\n",
    "\n",
    "\n",
    "guess_values = [1, t_0 - 10]\n",
    "q_1, t_2 = root(expressions, guess_values).x\n",
    "q_2 = (t_0 - t_1) * lambda_s / delta_s\n",
    "# 在薄膜与基板的接触面处能量守恒\n",
    "q = q_1 + q_2\n",
    "print(f'q = {q:.2f} W/m^2')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-13"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "相对误差大小为28.48%\n"
     ]
    }
   ],
   "source": [
    "Delta = 0.1e-3\n",
    "t_1 = 180\n",
    "t_2 = 30\n",
    "fluid = 'Air'\n",
    "Phi = 58.2\n",
    "d = 120e-3\n",
    "\n",
    "p = sc.atm\n",
    "T_1 = sc.convert_temperature(t_1, 'C', 'K')\n",
    "T_2 = sc.convert_temperature(t_2, 'C', 'K')\n",
    "\n",
    "lambda_1 = psi('L', 'T', T_1, 'P', p, fluid)\n",
    "lambda_2 = psi('L', 'T', T_2, 'P', p, fluid)\n",
    "\n",
    "area = np.pi * d**2 / 4    # 纵向单位长度的传热面积\n",
    "q = Phi / area\n",
    "R_total = (t_1 - t_2) / q\n",
    "R_1 = Delta / lambda_1\n",
    "R_2 = Delta / lambda_2\n",
    "R = R_total - R_1 - R_2\n",
    "\n",
    "error = R_total - R\n",
    "error_ratio = error / R\n",
    "\n",
    "print(f'相对误差大小为{error_ratio:.2%}')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-14"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "delta = 0.107 m\n"
     ]
    }
   ],
   "source": [
    "d = 100e-3\n",
    "density = 20\n",
    "t_o = 400\n",
    "t_insu = 50\n",
    "Phi = 163\n",
    "material = '超细玻璃棉毡、管'\n",
    "lambda_insu = get_lambda_(material, (t_o + t_insu) / 2)\n",
    "R_insu = (t_o - t_insu) / Phi\n",
    "lambda_insu = get_lambda_(material, (t_o + t_insu) / 2)\n",
    "\n",
    "guess_value = d\n",
    "delta = root(lambda delta: R_insu - cylindrical_wall_R(d/2, d/2 + delta, lambda_insu, 1), guess_value).x[0]\n",
    "print(f'delta = {delta:.3f} m')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-15"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "矿渣棉与煤灰泡沫砖交界界面处的温度不超过允许值\n",
      "增加煤灰泡沫砖的厚度会导致热损失下降\n"
     ]
    }
   ],
   "source": [
    "d_1 = 50e-3\n",
    "delta_1 = 40e-3\n",
    "lambda_1 = 0.11\n",
    "delta_2 = 45e-3\n",
    "lambda_2 = 0.12\n",
    "t_o = 50\n",
    "t_i = 400\n",
    "material1 = '矿渣棉'\n",
    "material2 = '粉煤灰泡沫砖'\n",
    "\n",
    "idx1 = ap4.material_list.index(material1)\n",
    "idx2 = ap4.material_list.index(material2)\n",
    "t_12_max = min(ap4.t_max[idx1], ap4.t_max[idx2])\n",
    "\n",
    "R_1 = cylindrical_wall_R(d_1/2, d_1/2+delta_1, lambda_1, 1)\n",
    "R_2 = cylindrical_wall_R(d_1/2+delta_1, d_1/2+delta_1+delta_2, lambda_2, 1)\n",
    "R = R_1 + R_2\n",
    "q = (t_i - t_o) / R\n",
    "t_12 = t_o + q * R_2\n",
    "if t_12 > t_12_max:\n",
    "    print('矿渣棉与煤灰泡沫砖交界界面处的温度超过允许值')\n",
    "else:\n",
    "    print('矿渣棉与煤灰泡沫砖交界界面处的温度不超过允许值')\n",
    "\n",
    "R_2, delta_2 = sp.symbols('R_2 delta_2')\n",
    "r2 = d_1/2 + delta_1 + delta_2\n",
    "r1 = d_1/2 + delta_1\n",
    "length = 1\n",
    "R_2 = sp.log(r2/r1) / (2*np.pi*lambda_2*length)\n",
    "dR_2_ddelta_2 = sp.diff(R_2, delta_2).subs(delta_2, 0.045)\n",
    "if dR_2_ddelta_2 > 0:\n",
    "    print('增加煤灰泡沫砖的厚度会导致热损失下降')\n",
    "else:\n",
    "    print('增加煤灰泡沫砖的厚度会导致热损失增加')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-16"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I = 232.42 A\n"
     ]
    }
   ],
   "source": [
    "d = 3e-3\n",
    "R_e = 2.22e-3 # Ohm per meter\n",
    "delta = 1e-3\n",
    "lambda_ = 0.15\n",
    "t_max = 65\n",
    "t_min = 0\n",
    "\n",
    "R = cylindrical_wall_R(d/2, d/2+delta, lambda_, 1)\n",
    "Phi = (t_max - t_min) / R\n",
    "I = np.sqrt(Phi / R_e)\n",
    "print(f'I = {I:.2f} A')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-17"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1) q = 12539.34 W\n",
      "(2) q = 5808.63 W\n",
      "(3) q = 10910.08 W\n"
     ]
    }
   ],
   "source": [
    "t_gas = 1000\n",
    "t_water = 200\n",
    "h_o = 100\n",
    "h_i = 5000\n",
    "delta = 6e-3\n",
    "lambda_ = 42\n",
    "d_o = 52e-3\n",
    "\n",
    "delta_dirty = [0, 1e-3, 2e-3]\n",
    "lamdba_dirty = [1, 0.08, 1]\n",
    "\n",
    "d_i = d_o - 2*delta\n",
    "d_o = d_o\n",
    "R_1 = 1 / (h_i * np.pi * d_i)\n",
    "R_2 = cylindrical_wall_R(d_i/2, d_o/2, lambda_, 1)\n",
    "\n",
    "for i, delta_d in enumerate(delta_dirty):\n",
    "    lambda_d = lamdba_dirty[i]\n",
    "    R_3 = cylindrical_wall_R(d_o/2, d_o/2 + delta_d, lambda_d, 1)\n",
    "    R_4 = 1 / (h_o * np.pi * (d_o + delta_d))\n",
    "    R_total = R_1 + R_2 + R_3 + R_4\n",
    "    q = (t_gas - t_water) / R_total\n",
    "    print(f'({i+1}) q = {q:.2f} W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-18"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "导热系数小的材料紧贴管壁保温效果更好\n"
     ]
    }
   ],
   "source": [
    "d = 100e-3\n",
    "lambda_1 = 0.05\n",
    "lambda_2 = 0.08\n",
    "delta_1 = delta_2 = 75e-3\n",
    "\n",
    "R_1i = cylindrical_wall_R(d/2, d/2+delta_1, lambda_1, 1)\n",
    "R_1o = cylindrical_wall_R(d/2+delta_1, d/2+delta_1+delta_2, lambda_1, 1)\n",
    "R_2i = cylindrical_wall_R(d/2, d/2+delta_2, lambda_2, 1)\n",
    "R_2o = cylindrical_wall_R(d/2+delta_2, d/2+delta_2+delta_1, lambda_2, 1)\n",
    "\n",
    "R_1 = R_1i + R_1o\n",
    "R_2 = R_2i + R_2o\n",
    "\n",
    "if R_1 > R_2:\n",
    "    print(f'导热系数小的材料紧贴管壁保温效果更好')\n",
    "else:\n",
    "    print(f'导热系数大的材料紧贴管壁保温效果更好')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-19"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "先案敷设材料B，再敷设材料A\n"
     ]
    }
   ],
   "source": [
    "d = 30e-3\n",
    "t = 100\n",
    "t_amb = 20\n",
    "Phi_l1 = 100   # W/m\n",
    "Phi_l2 = 50    # W/m, aim\n",
    "\n",
    "area = np.pi * d\n",
    "\n",
    "h = Phi_l1 / ((t - t_amb) * area)\n",
    "\n",
    "lambda_A = 0.5\n",
    "v_l_A = 3.14e-3\n",
    "lambda_B = 0.1\n",
    "v_l_B = 4.0e-3\n",
    "\n",
    "\n",
    "def find_outer_radius(r_in, area):\n",
    "    r_out = np.sqrt(area / np.pi + r_in**2)\n",
    "    return r_out\n",
    "\n",
    "\n",
    "# 第一种方案，先敷设材料A，再敷设材料B\n",
    "r_1 = find_outer_radius(d/2, v_l_A)\n",
    "r_2 = find_outer_radius(r_1, v_l_B)\n",
    "\n",
    "R_1 = cylindrical_wall_R(d/2, r_1, lambda_A, 1)\n",
    "R_2 = cylindrical_wall_R(r_1, r_2, lambda_B, 1)\n",
    "area_out = np.pi * r_2 * 2\n",
    "R_3 = 1 / (h * area_out)\n",
    "\n",
    "R = R_1 + R_2 + R_3\n",
    "Phi_l_A = (t - t_amb) / R\n",
    "\n",
    "# 第二种方案，先敷设材料B，再敷设材料A\n",
    "r_1 = find_outer_radius(d/2, v_l_B)\n",
    "r_2 = find_outer_radius(r_1, v_l_A)\n",
    "\n",
    "R_1 = cylindrical_wall_R(d/2, r_1, lambda_B, 1)\n",
    "R_2 = cylindrical_wall_R(r_1, r_2, lambda_A, 1)\n",
    "area_out = np.pi * r_2 * 2\n",
    "R_3 = 1 / (h * area_out)\n",
    "\n",
    "R = R_1 + R_2 + R_3\n",
    "Phi_l_B = (t - t_amb) / R\n",
    "\n",
    "if Phi_l_A < Phi_l2 and Phi_l_B < Phi_l2:\n",
    "    print('两种方案都可行')\n",
    "elif Phi_l_A < Phi_l2:\n",
    "    print('先案敷设材料A，再敷设材料B')\n",
    "elif Phi_l_B < Phi_l2:\n",
    "    print('先案敷设材料B，再敷设材料A')\n",
    "else:\n",
    "    print('两种方案都不可行')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-20"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1) t(x) = (l*t_1 - t_1*x + t_2*x)/l\n",
      "(2) t(x) = (t_1*exp(4*l*sqrt(h/(d*lambda_))) - t_1*exp(4*x*sqrt(h/(d*lambda_))) - t_2*exp(2*l*sqrt(h/(d*lambda_))) + t_2*exp(2*sqrt(h/(d*lambda_))*(l + 2*x)) - t_f*exp(4*l*sqrt(h/(d*lambda_))) + t_f*exp(2*l*sqrt(h/(d*lambda_))) + t_f*exp(4*x*sqrt(h/(d*lambda_))) - t_f*exp(2*x*sqrt(h/(d*lambda_))) - t_f*exp(2*sqrt(h/(d*lambda_))*(l + 2*x)) + t_f*exp(2*sqrt(h/(d*lambda_))*(2*l + x)))*exp(-2*x*sqrt(h/(d*lambda_)))/(exp(4*l*sqrt(h/(d*lambda_))) - 1)\n"
     ]
    }
   ],
   "source": [
    "d, l, t_1, t_2, t_f, lambda_, h = symbols('d, l, t_1, t_2, t_f, lambda_, h')\n",
    "x = symbols('x')\n",
    "t = Function('t')\n",
    "\n",
    "eq1 = diff(t(x), x, 2)\n",
    "result = dsolve(eq1, ics={t(0): t_1, t(l): t_2})\n",
    "t = solve(result, t(x))[0]\n",
    "print(f'(1) t(x) = {t}')\n",
    "\n",
    "t = Function('t')\n",
    "eq2 = diff(t(x), x, 2) - 4 * h / (d * lambda_) * (t(x) - t_f)\n",
    "result = dsolve(eq2, ics={t(0): t_1, t(l): t_2})\n",
    "t = solve(result, t(x))[0]\n",
    "print(f'(2) t(x) = {t}')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-21"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t(x) = (t_1*exp(4*l*sqrt(h/(d*lambda_))) - t_1*exp(4*x*sqrt(h/(d*lambda_))) - t_2*exp(2*l*sqrt(h/(d*lambda_))) + t_2*exp(2*sqrt(h/(d*lambda_))*(l + 2*x)) - t_f*exp(4*l*sqrt(h/(d*lambda_))) + t_f*exp(2*l*sqrt(h/(d*lambda_))) + t_f*exp(4*x*sqrt(h/(d*lambda_))) - t_f*exp(2*x*sqrt(h/(d*lambda_))) - t_f*exp(2*sqrt(h/(d*lambda_))*(l + 2*x)) + t_f*exp(2*sqrt(h/(d*lambda_))*(2*l + x)))*exp(-2*x*sqrt(h/(d*lambda_)))/(exp(4*l*sqrt(h/(d*lambda_))) - 1)\n",
      "钢柱体单位时间从两个热源分别获得的热量为:19.47 W 和2.01 W，总计21.48 W。\n"
     ]
    }
   ],
   "source": [
    "d, l, t_1, t_2, t_f, lambda_, h = symbols('d, l, t_1, t_2, t_f, lambda_, h')\n",
    "x = symbols('x')\n",
    "\n",
    "# eq1 = diff(t(x), x, 2)\n",
    "# result = dsolve(eq1, ics={t(0): t_1, t(l): t_2})\n",
    "# t = solve(result, t(x))[0]\n",
    "# print(f'(1) t(x) = {t}')\n",
    "\n",
    "t = Function('t')\n",
    "eq2 = diff(t(x), x, 2) - 4 * h / (d * lambda_) * (t(x) - t_f)\n",
    "result = dsolve(eq2, ics={t(0): t_1, t(l): t_2})\n",
    "t = solve(result, t(x))[0]\n",
    "print(f't(x) = {t}')\n",
    "\n",
    "d = 20e-3\n",
    "l = 300e-3\n",
    "t_1 = 250\n",
    "t_2 = 60\n",
    "t_f = 30\n",
    "\n",
    "h = 10\n",
    "lambda_ = 40\n",
    "\n",
    "area = np.pi * d ** 2 / 4\n",
    "phi = - lambda_ * diff(t, x, 1) * area\n",
    "phi_1 = phi.subs([(symbols('d'), d), (symbols('l'), l), (symbols('t_1'), t_1), (symbols('t_2'), t_2),\n",
    "                  (symbols('t_f'), t_f), (symbols('h'), h), (symbols('lambda_'), lambda_), (symbols('x'), 0)])\n",
    "phi_2 = phi.subs([(symbols('d'), 20e-3), (symbols('l'), l), (symbols('t_1'), t_1), (symbols('t_2'), t_2),\n",
    "                  (symbols('t_f'), 30), (symbols('h'), h), (symbols('lambda_'), lambda_), (symbols('x'), l)])\n",
    "phi_total = phi_1 + phi_2\n",
    "print(f'钢柱体单位时间从两个热源分别获得的热量为:{phi_1:.2f} W 和{phi_2:.2f} W，总计{phi_total:.2f} W。')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-22"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "上述条件下，液氮每天的蒸发量为：0.19 kg\n"
     ]
    }
   ],
   "source": [
    "d = 300e-3\n",
    "delta = 30e-3\n",
    "lambda_ = 1.8e-4\n",
    "t = -195.6\n",
    "t_a = 25\n",
    "gamma = 199.6e3\n",
    "time = 1*sc.day\n",
    "\n",
    "R = spherical_wall_R(d/2, d/2+delta, lambda_)\n",
    "Phi = (t_a - t) / R\n",
    "m = Phi * time / gamma\n",
    "print(f'上述条件下，液氮每天的蒸发量为：{m:.2f} kg')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-23"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 个球罐所必须配备的制冷设备的容量为：3547.00 W\n"
     ]
    }
   ],
   "source": [
    "d = 2\n",
    "delta = 1e-2\n",
    "t = -40\n",
    "delta_insu = 30e-2\n",
    "lambda_ = 0.08\n",
    "t_a = 40\n",
    "h_insu = 30\n",
    "number = 10\n",
    "\n",
    "r_1 = d/2 + delta   # 认为绝缘层内壁温度等于罐内温度，只分析绝缘层的传热\n",
    "r_2 = d/2 + delta + delta_insu\n",
    "R = spherical_wall_R(r_1, r_2, lambda_)\n",
    "Phi = number * (t_a - t) / R\n",
    "print(f'10 个球罐所必须配备的制冷设备的容量为：{Phi:.2f} W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-24"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lambda = 1.20 W/m-K\n"
     ]
    }
   ],
   "source": [
    "d_i = 0.15\n",
    "d_o = 0.25\n",
    "t_i = 200\n",
    "t_o = 40\n",
    "P = 56.5\n",
    "\n",
    "R = P / (t_i - t_o)\n",
    "guess_value = 1\n",
    "lambda_ = root(lambda lambda_: R - spherical_wall_R(d_i/2, d_o/2, lambda_), guess_value).x[0]\n",
    "print(f'lambda = {lambda_:.2f} W/m-K')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "outputs": [],
   "source": [
    "## 习题02-25"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "球罐的内表面温度为：53.19 C\n",
      "球罐的外表面温度为：30.79 C\n"
     ]
    }
   ],
   "source": [
    "d_i = 0.5\n",
    "d_o = 0.6\n",
    "Phi_dot = 1e5   # W/m^3\n",
    "h = 1000        # W/m^2-K\n",
    "t_f = 25\n",
    "material = '铬镍钢'\n",
    "\n",
    "temperature_list = [-100, 0, 100, 200, 300, 400, 600, 800]\n",
    "lambda_list = [12.2, 14.7, 16.6, 18.0, 19.4, 20.8, 23.5, 26.3]\n",
    "\n",
    "\n",
    "def lambda_(t):\n",
    "    return np.interp(t, temperature_list, lambda_list)\n",
    "\n",
    "\n",
    "V = 4/3 * np.pi * (d_i/2)**3\n",
    "Phi = Phi_dot * V\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    t_i, t_o = p\n",
    "    xpr1 = Phi - (t_i - t_o) / spherical_wall_R(d_i/2, d_o/2, lambda_((t_i+t_o)/2))\n",
    "    xpr2 = Phi - h * np.pi * d_o**2 * (t_o - t_f)\n",
    "    return [xpr1, xpr2]\n",
    "\n",
    "\n",
    "guess_values = [t_f + 20, t_f + 10]\n",
    "t_i, t_o = root(expressions, guess_values).x\n",
    "print(f'球罐的内表面温度为：{t_i:.2f} C')\n",
    "print(f'球罐的外表面温度为：{t_o:.2f} C')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-26"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "所需的电加热器的功率为：5198.15 W\n"
     ]
    }
   ],
   "source": [
    "delta = 20e-2\n",
    "lambda_ = 1.5\n",
    "\n",
    "T_i = 400\n",
    "t_f = 25\n",
    "h = 10\n",
    "r_o = 0.5\n",
    "l = 2.0\n",
    "\n",
    "t_i = sc.convert_temperature(T_i, 'K', 'C')\n",
    "\n",
    "# 可将导热热阻开成球壳热阻和圆柱壳热阻的并联\n",
    "# 圆柱壳热阻\n",
    "R_c = cylindrical_wall_R(r_o - delta, r_o, lambda_, l)\n",
    "# 球壳热阻\n",
    "R_s = spherical_wall_R(r_o - delta, r_o, lambda_)\n",
    "\n",
    "R = 1 / (1 / R_c + 1 / R_s)\n",
    "Phi = (t_i - t_f) / R\n",
    "print(f'所需的电加热器的功率为：{Phi:.2f} W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-27"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Phi = 0.036 W\n",
      "Phi_p = 0.046 W\n"
     ]
    }
   ],
   "source": [
    "r_1 = 10e-3\n",
    "r_2 = 12.5e-3\n",
    "r_3 = 16.3e-3\n",
    "t_fi = 37\n",
    "t_fo = 20\n",
    "h_i = 12\n",
    "h_o = 6\n",
    "lambda_1 = 0.35\n",
    "lambda_2 = 0.8\n",
    "ratio = 1/3\n",
    "\n",
    "area_1 = ratio * 4 * np.pi * r_1**2\n",
    "area_2 = ratio * 4 * np.pi * r_2**2\n",
    "area_3 = ratio * 4 * np.pi * r_3**2\n",
    "\n",
    "R_i = 1 / (h_i * area_1)\n",
    "# 由于球面可看做各部分热阻并联组成，所以热阻等于球壳的热阻除以ratio\n",
    "R_1 = spherical_wall_R(r_1, r_2, lambda_1) / ratio\n",
    "R_2 = spherical_wall_R(r_2, r_3, lambda_2) / ratio\n",
    "R_o = 1 / (h_o * area_2)\n",
    "R_o_p = 1 / (h_o * area_3)\n",
    "\n",
    "R_total = R_i + R_1 + R_o\n",
    "R_total_p = R_i + R_1 + R_2 + R_o_p\n",
    "\n",
    "Phi = (t_fi - t_fo) / R_total\n",
    "Phi_p = (t_fi - t_fo) / R_total_p\n",
    "\n",
    "print(f'Phi = {Phi:.3f} W')\n",
    "print(f'Phi_p = {Phi_p:.3f} W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-28"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "结冰的冰软木保温层的厚度为：0.301 m\n"
     ]
    }
   ],
   "source": [
    "d = 1.22\n",
    "delta = 0.45\n",
    "lambda_ = 0.043\n",
    "t_i = -62.2\n",
    "t_o = 18\n",
    "h_i = 1050\n",
    "h_o = 21\n",
    "\n",
    "#   做出以下假设：\n",
    "#   1. 结冰的软木保温层由含水层和干软木保温层组成\n",
    "#   2. 含水层的导热系数取冰的导热系数\n",
    "#   3. 干软木保温层的导热系数取软木的导热系数\n",
    "#   4. 含水（冰）软木保温层的厚度与之前相同\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    delta_ice = p[0]\n",
    "    r_1 = d / 2\n",
    "    r_2 = d / 2 + delta_ice # delta_ice <= delta\n",
    "    r_3 = d / 2 + delta\n",
    "\n",
    "    area_i = 4 * np.pi * r_1**2\n",
    "    area_o = 4 * np.pi * r_3**2\n",
    "\n",
    "    R_i = 1 / (h_i * area_i)\n",
    "    R_ice = spherical_wall_R(r_1, r_2, lambda_)\n",
    "    R_wood = spherical_wall_R(r_2, r_3, lambda_)\n",
    "    R_o = 1 / (h_o * area_o)\n",
    "\n",
    "    R = R_i + R_ice + R_wood + R_o\n",
    "    Phi = (t_o - t_i) / R\n",
    "    t_2 = 0\n",
    "    xpr = Phi - (t_o - t_2)/(R_wood + R_o)\n",
    "    return xpr\n",
    "\n",
    "\n",
    "guess_value = delta/2\n",
    "delta_ice = root(expressions, guess_value).x[0]\n",
    "print(f'结冰的冰软木保温层的厚度为：{delta_ice:.3f} m')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-29"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t = 51.53 C\n"
     ]
    }
   ],
   "source": [
    "r = 0.1e-3\n",
    "t_oo = 25\n",
    "lambda_ = 120\n",
    "Phi = 4\n",
    "\n",
    "R = spherical_wall_R(r, np.inf, lambda_)\n",
    "guess_value = t_oo + 1\n",
    "t = root(lambda t: Phi - (t - t_oo) / R, guess_value).x[0]\n",
    "print(f't = {t:.2f} C')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-30"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Phi = 1445.30 W\n"
     ]
    }
   ],
   "source": [
    "height = 30e-2\n",
    "d_1 = 8.2e-2\n",
    "d_2 = 13e-2\n",
    "t_1 = 20\n",
    "t_2 = 520\n",
    "lambda_ = 100\n",
    "\n",
    "A = Function('A')\n",
    "t = symbols('t')\n",
    "Phi = symbols('Phi')\n",
    "y = symbols('y')\n",
    "\n",
    "eq1 = diff(A(y), y, 2)\n",
    "A = solve(dsolve(eq1, ics={A(0): pi/4*d_2**2, A(height): pi/4*d_1**2}), A(y))[0]\n",
    "# Phi = -lambda_*A*diff(t, y)\n",
    "result1 = integrate(-1/(lambda_ * A), (y, 0, height))\n",
    "result2 = integrate(1, (t, t_2, t_1))\n",
    "Phi = result2 / result1\n",
    "print(f'Phi = {Phi:.2f} W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-31"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "第3种导热问题热流量排第1\n",
      "第1种导热问题热流量排第2\n",
      "第2种导热问题热流量排第3\n"
     ]
    }
   ],
   "source": [
    "a_list = np.array([0.506, 0.08, 20.24])\n",
    "n_list = np.array([0.5, 0.0, 1.5])\n",
    "x_1 = 25e-3\n",
    "x_2 = 125e-3\n",
    "\n",
    "number = len(a_list)\n",
    "result = np.zeros(number)\n",
    "for i in range(number):\n",
    "    a = a_list[i]\n",
    "    n = n_list[i]\n",
    "\n",
    "    def func(x):\n",
    "        d = a*x**n\n",
    "        area = np.pi / 4 * d**2\n",
    "        return 1/area\n",
    "    result[i] = 1 / quad(func, x_1, x_2)[0]\n",
    "\n",
    "arg = np.argsort(-result)\n",
    "for i in range(number):\n",
    "    print(f'第{arg[i]+1}种导热问题热流量排第{i+1}')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-32"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lambda = 5.06\n",
      "lambda_0 = 5.05555555555556*(2.0*t**2 - 8825.0)/((t - 85.0)*(t - 40.0))\n",
      "b = -2.0*(2.0*t - 125.0)/(2.0*t**2 - 8825.0)\n"
     ]
    }
   ],
   "source": [
    "delta = 25e-3\n",
    "t_1 = 40\n",
    "t_2 = 85\n",
    "phi = 1.82e3\n",
    "area = 0.2\n",
    "\n",
    "lambda_ = phi * delta / (area * (t_2 - t_1))\n",
    "print(f'lambda = {lambda_:.2f}')\n",
    "\n",
    "# 如果已知平板中间的温度t，则可以得到lmbda的公式\n",
    "\n",
    "lambda_1 = phi * delta / (area * (t - t_1))\n",
    "lambda_2 = phi * delta / (area * (t_2 - t))\n",
    "\n",
    "lambda_0, b = symbols('lambda_0 b')\n",
    "\n",
    "result = solve([lambda_1 - lambda_0 * (1 + b*(t+t_1)/2), lambda_2 - lambda_0 * (1 + b*(t+t_2)/2)], lambda_0, b)\n",
    "print(f'lambda_0 = {result[0][0]}')\n",
    "print(f'b = {result[0][1]}')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-36"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lambda_0 = 0.917 W/m-K\n",
      "b = -9.091e-03 K^-1\n"
     ]
    }
   ],
   "source": [
    "q = 1000\n",
    "delta = 20e-3\n",
    "x_list = np.array([0, 10e-3, 20e-3])\n",
    "t_list = np.array([100, 60, 40])\n",
    "\n",
    "\n",
    "# lambda_ = lambda_0 * (1 + b*t)\n",
    "\n",
    "\n",
    "def lambda_(delta, delta_t):\n",
    "    return -q * delta / delta_t\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    lambda_0, b = p\n",
    "\n",
    "    def lambda_b(t):\n",
    "        return lambda_0 * (1 + b * t)\n",
    "\n",
    "    delta_1 = x_list[1] - x_list[0]\n",
    "    delta_2 = x_list[2] - x_list[1]\n",
    "    delta_t_1 = t_list[1] - t_list[0]\n",
    "    delta_t_2 = t_list[2] - t_list[1]\n",
    "    t_average_1 = (t_list[1] + t_list[0]) / 2\n",
    "    t_average_2 = (t_list[2] + t_list[1]) / 2\n",
    "    xpr1 = lambda_(delta_1, delta_t_1) - lambda_b(t_average_1)\n",
    "    xpr2 = lambda_(delta_2, delta_t_2) - lambda_b(t_average_2)\n",
    "    return xpr1, xpr2\n",
    "\n",
    "\n",
    "guess_values = [1, 1]\n",
    "lambda_0, b = root(expressions, guess_values).x\n",
    "print(f'lambda_0 = {lambda_0:.3f} W/m-K')\n",
    "print(f'b = {b:.3e} K^-1')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-41"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "最大热功率为：77875.46 W\n"
     ]
    }
   ],
   "source": [
    "t_max = 1600\n",
    "t_water = 110\n",
    "h = 12000\n",
    "R_c = 2.2e-4    # m^2-K/W\n",
    "r_1 = 6.1e-3\n",
    "r_2 = 6.5e-3\n",
    "lambda_ = 7.9\n",
    "lambda_alloy = 14.2\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    phi = p\n",
    "    q_dot = phi / (np.pi * r_1**2)\n",
    "    t_w = t_max - cylinder_Delta_T_with_heat_source(r_1, lambda_, q_dot)\n",
    "    R_w = cylindrical_wall_R(r_1, r_2, lambda_, 1)\n",
    "    R_conv = 1 / (h * 2*np.pi*r_2)\n",
    "    R = R_c / (2*np.pi*r_1) + R_w + R_conv  # R_c的单位表明其需要除以接触面积\n",
    "    xpr = phi - (t_w - t_water) / R\n",
    "    return xpr\n",
    "\n",
    "\n",
    "guess_values = 1e6\n",
    "phi = root(expressions, guess_values).x[0]\n",
    "print(f'最大热功率为：{phi:.2f} W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-46"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t_max = 117.50 C\n"
     ]
    }
   ],
   "source": [
    "# 绝热表面可以看做对称面\n",
    "delta = 7e-2\n",
    "t_f = 30\n",
    "q_dot = 0.3e6\n",
    "h = 450\n",
    "lambda_ = 18\n",
    "\n",
    "Delta_T = plate_wall_Delta_T_with_heat_source(2*delta, lambda_, h, q_dot)\n",
    "t_max = t_f + Delta_T\n",
    "print(f't_max = {t_max:.2f} C')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-50"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1) eta = 92.43%\n",
      "(2) eta = 60.33%\n"
     ]
    }
   ],
   "source": [
    "lambda_ = np.array([208, 41.5])\n",
    "h = np.array([284, 511])\n",
    "H = np.array([15.24e-3, 15.24e-3])\n",
    "delta = np.array([2.54e-3, 2.54e-3])\n",
    "\n",
    "perimeter = 2\n",
    "A_c = delta\n",
    "efficiency = fin_tip_efficiency(H, perimeter, A_c, lambda_, h)\n",
    "# 也可以用fin_tip_efficiency2函数求解\n",
    "# A_L = delta * H\n",
    "# efficiency2 = fin_tip_efficiency2(H, A_L, lambda_, h)\n",
    "np.set_printoptions(formatter={'float': '{:.2%}'.format})\n",
    "for i, eta in enumerate(efficiency):\n",
    "    print(f'({i+1}) eta = {eta:.2%}')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-51"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Phi2 = 39.74 W\n",
      "Phi2 = 65.12 W\n"
     ]
    }
   ],
   "source": [
    "material = 'Al'\n",
    "t_w = 260\n",
    "d = 25e-3\n",
    "H = 150e-3\n",
    "t_1 = 16\n",
    "h = 15\n",
    "\n",
    "perimeter = np.pi * d\n",
    "A_c = np.pi * d**2 / 4\n",
    "lambda_ = 208   # 采用上一题给出Al的导热系数\n",
    "\n",
    "R = fin_tip_R(H, perimeter, A_c, lambda_, h)\n",
    "Phi = (t_w - t_1) / R\n",
    "print(f'Phi2 = {Phi:.2f} W')\n",
    "\n",
    "H_2 = 2*H\n",
    "R_2 = fin_tip_R(H_2, perimeter, A_c, lambda_, h)\n",
    "Phi_2 = (t_w - t_1) / R_2\n",
    "print(f'Phi2 = {Phi_2:.2f} W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-52"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mH = 0.4474091316149847\n",
      "总散热量为：4541.22 W\n"
     ]
    }
   ],
   "source": [
    "material = 'Al'\n",
    "d = 25e-3\n",
    "s = 9.5e-3\n",
    "H = 12.5e-3\n",
    "delta = 0.8e-3\n",
    "t_w = 200\n",
    "t_f = 90\n",
    "h = 110\n",
    "lambda_ = 208   # 参考题2-50取值\n",
    "length = 1\n",
    "\n",
    "H_p = H + delta/2\n",
    "A_L = H_p * delta\n",
    "fin_efficiency = fin_tip_efficiency2(H, A_L, lambda_, h)\n",
    "fin_area = 2 * np.pi * ((d/2+H_p)**2 - (d/2)**2)\n",
    "Phi_fin = fin_efficiency * h * fin_area * (t_w - t_f)\n",
    "base_area = np.pi * d * (s - delta)\n",
    "Phi_base = h * base_area * (t_w - t_f)\n",
    "number = length / s\n",
    "Phi = number * (Phi_fin + Phi_base)\n",
    "\n",
    "print(f'总散热量为：{Phi:.2f} W')\n",
    "\n",
    "# 也可以采用overall_fin_surface_efficiency函数求解\n",
    "# eta = overall_fin_surface_efficiency(fin_area, base_area, fin_efficiency)\n",
    "# Phi = eta * h * (fin_area + base_area) * (t_w - t_f) * number\n",
    "# print(f'总散热量为：{Phi:.2f} W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-53"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "套管应有的长度为：0.123 m\n"
     ]
    }
   ],
   "source": [
    "D = 127e-3\n",
    "d = 15e-3\n",
    "delta = 0.9e-3\n",
    "lambda_ = 49.1\n",
    "h = 105\n",
    "Delta_T_ratio = 0.6e-2\n",
    "\n",
    "perimeter = np.pi * d\n",
    "A_c = np.pi * (d/2 + delta)**2 - np.pi * (d/2)**2\n",
    "\n",
    "guess_value = 0.1\n",
    "H = root(lambda H: Delta_T_ratio - fin_tip_Delta_T_ratio(H, perimeter, A_c, lambda_, h), guess_value).x[0]\n",
    "print(f'套管应有的长度为：{H:.3f} m')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-54"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1) Delta_T = 31.04 C\n",
      "(2) Delta_T = 6.93 C\n",
      "两温度计读数相差24.11 C\n"
     ]
    }
   ],
   "source": [
    "d = 10e-3\n",
    "delta = 1.0e-3\n",
    "H = 120e-3\n",
    "h = 25\n",
    "t_0 = 25\n",
    "t_f = 70\n",
    "lambda_ = [390, 50]\n",
    "\n",
    "perimeter = np.pi * d\n",
    "A_c = np.pi * (d/2 + delta)**2 - np.pi * (d/2)**2\n",
    "Delta_T = np.zeros(2)\n",
    "\n",
    "for i, lambda_i in enumerate(lambda_):\n",
    "    Delta_T_ratio = fin_tip_Delta_T_ratio(H, perimeter, A_c, lambda_i, h)\n",
    "    Delta_T[i] = Delta_T_ratio * (t_f - t_0)\n",
    "    print(f'({i+1}) Delta_T = {Delta_T[i]:.2f} C')\n",
    "\n",
    "delta_T12 = Delta_T[0] - Delta_T[1]\n",
    "print(f'两温度计读数相差{abs(delta_T12):.2f} C')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-55"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "该柱体中间截面上的平均温度为392.71 C\n",
      "该柱体的最高温度温度为548.98 C\n",
      "冷却介质带走的热量为为65.74 W\n"
     ]
    }
   ],
   "source": [
    "H = 9e-2\n",
    "perimeter = 7.6e-2\n",
    "A_c = 1.95e-4\n",
    "t_0 = 305\n",
    "t_f = 815\n",
    "h = 28\n",
    "lambda_ = 55\n",
    "\n",
    "y = H/2\n",
    "\n",
    "Delta_T_ratio_y = fin_tip_Delta_T_ratio(y, perimeter, A_c, lambda_, h)\n",
    "t_y = t_f - Delta_T_ratio_y*(t_f - t_0)\n",
    "print(f'该柱体中间截面上的平均温度为{t_y:.2f} C')\n",
    "\n",
    "Delta_T_ratio_H = fin_tip_Delta_T_ratio(H, perimeter, A_c, lambda_, h)\n",
    "t_H = t_f - Delta_T_ratio_H*(t_f - t_0)\n",
    "print(f'该柱体的最高温度温度为{t_H:.2f} C')\n",
    "\n",
    "efficiency = fin_tip_efficiency(H, perimeter, A_c, lambda_, h)\n",
    "Phi = efficiency * h * perimeter * H * (t_f - t_0)\n",
    "print(f'冷却介质带走的热量为为{Phi:.2f} W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-56"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "手柄的散热量为2.88 W\n",
      "手柄的最低温度为24.85 C\n",
      "手柄的导热系数越大，散热量越大，手柄每一截面处的温度也越均匀。\n"
     ]
    }
   ],
   "source": [
    "d = 25e-3\n",
    "theta = np.pi\n",
    "R = 75e-3\n",
    "t_w = 80\n",
    "t_f = 20\n",
    "h = 10\n",
    "lambda_ = 1.5\n",
    "\n",
    "R_p = R - d/2\n",
    "H = theta * R_p / 2\n",
    "r = d/2\n",
    "A_c = np.pi * r**2\n",
    "perimeter = np.pi * d\n",
    "\n",
    "R = fin_tip_R(H, perimeter, A_c, lambda_, h)\n",
    "Phi = 2 * (t_w - t_f) / R\n",
    "print(f'手柄的散热量为{Phi:.2f} W')\n",
    "\n",
    "eta = fin_tip_efficiency(H, perimeter, A_c, lambda_, h)\n",
    "Delta_T_ratio = fin_tip_Delta_T_ratio(H, perimeter, A_c, lambda_, h)\n",
    "t_min = t_f + Delta_T_ratio * (t_w - t_f)\n",
    "print(f'手柄的最低温度为{t_min:.2f} C')\n",
    "print(f'手柄的导热系数越大，散热量越大，手柄每一截面处的温度也越均匀。')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-57"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mH = 0.3606193232430791\n",
      "增加了肋片后汽缸散热量是原来的8.81倍\n"
     ]
    }
   ],
   "source": [
    "Diameter = 60e-3\n",
    "Height = 170e-3\n",
    "lambda_ = 180\n",
    "number = 10\n",
    "delta = 3e-3\n",
    "H = 25e-3\n",
    "h = 50\n",
    "t_f = 28\n",
    "t_w = 220\n",
    "\n",
    "Phi_0 = h * np.pi * Diameter * Height * (t_w - t_f)\n",
    "H_p = H + delta/2\n",
    "A_L = H_p * delta\n",
    "efficiency = fin_tip_efficiency2(H_p, A_L, lambda_, h)\n",
    "fin_tip_area = 2 * np.pi * ((Diameter/2 + H)**2 - (Diameter/2)**2)\n",
    "Phi_tip = number * efficiency * h * 2 * fin_tip_area *(t_w - t_f)\n",
    "Phi_base = h * np.pi * Diameter * (Height - number * delta) * (t_w - t_f)\n",
    "Phi_1 = Phi_tip + Phi_base\n",
    "times = Phi_1/ Phi_0\n",
    "print(f'增加了肋片后汽缸散热量是原来的{times:.2f}倍')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-58"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEECAYAAAAoDUMLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAApHElEQVR4nO3deXiU9dXG8e9JAmEJa4CA7IgYRIqQsNSKLIJibbWKtlaLW4UWaxdrqbbVt7a+ttbaxVdbq4iiVo3Yqi2KS1WiKAgEEEH2XRAIa0gCWee8f8xQYwzbJJPJzNyf65qLeeaZ38w5JJl7nt3cHRERkaRoFyAiIg2DAkFERAAFgoiIhCgQREQEUCCIiEhISrQLqI127dp5jx49whpbXFxM8+bN67agBk49J4ZE6znR+oXa97xo0aLd7t6++uMxHQg9evQgLy8vrLG5ubmMHDmybgtq4NRzYki0nhOtX6h9z2a2uabHtcpIREQABYKIiIQoEEREBFAgiIhIiAJBREQABYKIiIQoEEREBIjx4xBERGKRu3OovJLCkgoKS8opKq3kYGkFB8sqOVheSWl5JSUVAcoqAlRUBqgIOJUBJ+BOwGHjpjLOGh4gJbluv9MrEERE6kAg4OwpLmNHQQk7D5Sws7CE/AOl7CkuZU9RGXuKyth3sIx9B8spOFRGeWX416Ix4PcBJyW57uoHBYKIyHFxd/YWl7F570E27ylmy55DbN13kK37DrF1/0F2FpRSVhn43Lg2zRqRnpZK2+aNObl9Gm2aN6JV08a0atqIFk1SaNEkhbTUFJqnptCscTJNGyXTpFEyqY2SSE1OJiXZSEk2ks1ITjLMjNzcXJo0quM0QIEgIvIZlQHn470HWZtfxNr8QtbtLGLD7mI27CriQEnFZ56b0TKVLm2aMbBrG07q35ROrZrQsVUTOrZsQoeWqbRLS6VRHa/WiSQFgogkrMKSclZuL2TFJwWs2H6A1TsKWb2zkJLyT7/pZ7RM5eT2aXx1wEn0ap9Gj/RmdE9vRpc2zSLyLT2aFAgikhAOlVWybFsBH27dz4dbC1i2rYCNu4v/Oz+9eWMyO7XgiiHdyezYglMy0ji5QxotmzSKYtX1S4EgInHH3dm67xCLt+wjb9M+Fm/Zx6odhVQGghtyT2rVhP5dWjF+UGf6ndSKfie1pH2LVMwsypVHlwJBRGKeu7N+VzHzN+5h/oa9LNi4lx0HSgBo3jiZAV1b890RvRjYtQ0DuramfYvUKFfcMCkQRCQmbd13kPfW7eaFpSVMee9NdhWWAtC+RSpDe7ZlcI+2ZHVvQ2bHFnW+v368UiCISEwoLq1g3vo9vL1mF++u2/3f9f+tUo0Rmel88eR0hvVKp0d6s4Rf9RMuBYKINFgbdxfz5sqdzF6dz8KN+yirDNCscTLDeqUzYVh3hp/Sjq0r8hg1amC0S40LCgQRaTAqA86izfv4z4odvLkynw2hpYBTOqRxzZd6MLJPe7J6tCG1yiG621ZqaaCuKBBEJKpKKyqZu24Pry7fwX9W7mRvcRmNko1hvdK5+swejM7sQNe2zaJdZkKIWCCY2RTgQqAlcD+wErg3NPttd7+12vNbAhuA5aGHXnD3+yJVn4hET2lFJXPW7OblZdt5Y8VOCksraJGawqjMDpzbL4MRfdrTIoH2/28oIhIIZjYYGA6MANKAW4DJwHh332RmuWaW5e6LqgwbBMxw9xsiUZOIRFdFZYD31u/h3x98wusf7aCwtIJWTRsx7vSOfLl/J87snf6ZVUFS/8w9/DPuHfFFze4AmgGZfBoIS9y9wszSgJnApe6+p8qYm4FLgTJgJ/BDd99ew2tPAiYBZGRkZOXk5IRVY1FREWlpaWGNjVXqOTE0pJ7dnQ0FAeZ+UsGCHRUUlkHTFMjOSGFwx2ROS08mJal22wAaUr/1pbY9jxo1apG7Z1d/PFKB8BDQBxgHdCcYAJnAUOAZIB/4WtUPfDP7KlDq7q+b2ZUElyYuOdr7ZGdne15eXlg15ubmMnLkyLDGxir1nBgaQs9b9x3khcXbeH7JNjbuLiY1JYkxfTO48IyTGHlq+zpdEmgI/da32vZsZjUGQqS2IZQAs9y9FFhjZkVAO3d/H+gZWoL4MTClypi3QuMAXgB+FaHaRCQCSsoreWX5dp7L28rc9cGF/2G92jJ5xMmM698xoc4JFKsiFQjzgGvN7A9ABsHVRi+a2QXuvh8o5NMP/8OmElySeAYYCyyIUG0iUoeWbysgZ+EW/vXBJxSWVNCtbTNuGtOHSwZ11t5BMSZSgTADyALmhqYnA62A18zsELAZ+I6ZtQf+6u6XAT8HHjWz7wBFwMQI1SYitVRcWsG/l37CMwu28OHWAlJTkvhy/058PbsrQ3u2JamW2wUkOiISCO4e4LOrgw57odp0CXBZaMwmYHQk6hGRurFmZyF/f38zzy/eRlFpBadmtOBXF/bjawM706qpVgnFOh2YJiJHVVEZ4PUVO5k+dxMLNu6lcUoSX+nfiSuHdWNQtzY6b1AcUSCISI32FZfxzMIt/H3eZj4pKKFLm6bcen4mX8/uStvmjaNdnkSAAkFEPmNdfhGPvreR5xdvpaQ8wJd6p/Ori05ndGYHkrVtIK4pEEQEd2f+xr1MfWcDb67Kp3FKEhef0ZnrzurJqR1bRLs8qScKBJEEVhlwXl2+g4feWc+HWwtIb96Ym8b04VvDupGepquKJRoFgkgCKimv5PnF23j4nfVs2nOQnu2ac9fFpzN+UBeaNNL5hBKVAkEkgRwsq+Dp+Vt4+J0N5BeWMqBLKx68chDn9uuo7QOiQBBJBIUl5TwxbzOPzNnAvoPlnHlyOn/6xhmceXK6dhuV/1IgiMSxgkPlTH9vE9Pe3cCBkgpGZ3bgxtG9GdStTbRLkwZIgSAShwpLynnsvU08MicYBGNPy+AHo0+hf5dW0S5NGjAFgkgcOVhWwUsbyvjh27MpOFTO2NMy+OE5p3B6ZwWBHJsCQSQOlJRX8tT8LTyYu47dReWMzuzATWP6aIlATogCQSSGVVQGeH7JNv78nzV8UlDCl3qn8930Yq6/eHC0S5MYpEAQiUHuzusrdvL711azLr+IAV1a8fvLBvCl3u3Izc2NdnkSoxQIIjFm0eZ9/HbWSvI276NX++b87VuDOK9fR+0+KrWmQBCJEZt2F/O7V1fxyvIdtG+Rym8u7s/Xs7uQkpwU7dIkTigQRBq4/QfL+L831/Hk+5tolJzEj8f24frhPWnWWH++Urf0GyXSQJVXBnjq/c386Y21FJaU843BXblpbB86tGgS7dIkTkUsEMxsCnAh0BK4H1gJ3Bua/ba731rt+U2AqUBX4BBwtbvnR6o+kYZs9up8/velFazfVcxZvdtx21f6ktmxZbTLkjgXkUAws8HAcGAEkAbcAkwGxrv7JjPLNbMsd19UZdi1wBZ3n2BmVwC/AH4YifpEGqqNu4u586UVvLUqn57tmvPIVdmc07eDNhhLvTB3r/sXNbsDaAZk8mkgLHH3CjNLA2YCl7r7nipjngEecPf3zKw18Ja7D6rhtScBkwAyMjKycnJywqqxqKiItLS0sMbGKvXccJVUOP9eX85rm8pplAQX9W7M2O4ppIRxBtJY6bmuJFq/UPueR40atcjds6s/HqlVRp2APsA4oDvBAMg0s2HAM0A+UP2irOnA/tD9otD057j7w8DDANnZ2T5y5MiwCszNzSXcsbFKPTc87s5LH27nrpdXsuNAOeMHdeGW80+t1XaCht5zXUu0fiFyPUcqEEqAWe5eCqwxsyKgnbu/D/QMLUH8GJhSZcw+4PC1+loDexCJY+vyC/mff33E3PV76HdSS/5y5SCyuusspBI9kQqEecC1ZvYHIIPgaqMXzewCd98PFBIMjarmAOcB7wMXhKZF4s6hskruf2stU+dsoGmjZO782ulcMaSbLlAjURepQJgBZAFzQ9OTgVbAa2Z2CNgMfMfM2gN/dffLgGnAY2aWSzAsro5QbSJRM3tVPre9uJxt+w8xflAXfvblTNrp2sXSQEQkENw9wGdXBx32QrXpEuCy0JhDwOWRqEck2nYeKOHXM1fw8rLt9O6QxrOThjG0V42byUSiRgemiURQIOA8tWAL97yyirLKAFPOO5WJw3vROEWnm5CGR4EgEiFrdxbys+eXkbd5H1/qnc5dX+tPj3bNo12WyBEpEETqWFlFgL+9vZ7731pL89QU7r1sAOMHddbBZdLgKRBE6tDSj/dzyz8/ZNWOQr464CR++dXTtNFYYoYCQaQOlJRX8uc31vLwO+tp3yKVqVdlM/a0jGiXJXJCFAgitbRo8z5++o+lrN9VzOWDu/LzC/rSskmjaJclcsIUCCJhKimv5E//WcPUORvo1KopT1w3hLP7tI92WSJhUyCIhOHDrfu5ecZS1uYX8c0h3fj5lzNpoaUCiXEKBJETUF4Z4P631vGX2eton5bK49cNYYSWCiROKBBEjtPanYXcNOMDlm87wCUDO/PLC/vRqqmWCiR+KBBEjiEQcKbP3cTdr64iLTWFv31rEONO7xTtskTqnAJB5Ch2HijhJ88tZc7a3YzO7MDd4/vrmsYStxQIIkfw6vLt3Pr8MkrLA9x1cfAU1TraWOKZAkGkmuLSCn49cwXP5n3MF7q04s/fOINe7RPrEo2SmBQIIlUs21rAD3KWsGlPMTeMPJmbxvahUbLOTCqJQYEgQnDD8bR3N3LPa6tIb57K09cP44sn63oFklgUCJLwdheVcvOMpby9ZhfnnpbBPZd+gdbNGke7LJF6p0CQhDZ33W5++OwHFBwq586L+vGtYd214VgSVsQCwcymABcCLYH7gfXAXUAFkA9MCF028/DzWwIbgOWhh15w9/siVZ8ktsqAc9+ba7n/rbX0atecJ64bQt9OLaNdlkhURSQQzGwwMBwYAaQBtxC8xvJId99uZr8DricYFIcNAma4+w2RqEnksPwDJfwgZwnvb9jL+EFd+PVF/WieqoVlkUj9FVwArAJe5NNA+Ku7bw/NN4JLClVlAQPN7G1gJ/DDKs8XqRPvrt3Nj55dQnFpJX+4bADjs7pEuySRBsPcve5f1OwhoA8wDugOzAQy3d3NbDzwM+Bsdz9YZcxXgVJ3f93MrgTGu/slNbz2JGASQEZGRlZOTk5YNRYVFZGWllj7lidyzwF3/rWunH+vL6dTmvG9M5rQOS0+dydNtJ9zovULte951KhRi9w9+3Mz3L3Ob8B9wJQq04uA9sCNwHtAhxrGNAeSQ/ebAeuO9T5ZWVkertmzZ4c9NlYlas97ikr9W4+8791veclvenaJF5eWR7usiEq0n3Oi9ete+56BPK/hMzVSX5HmAWPMLMnMOhFcbTQJGAOc4+75NYyZCnw9dH8ssCBCtUkCWbe/kgv+bw7zN+7lt5f05w+XDaBZY20vEKlJpP4yZhDcJjA3NH078BTBJYVXQ7v1PQv8g+C2hcuAnwOPmtl3gCJgYoRqkwTg7jz5/mZ+O7+Ek9o05fnJZ3J651bRLkukQYtIILh7gOBeRVXNOMLTLwuN2QSMjkQ9klgOlVXy8xeW8cKSbQxon8wTk4fTqpmuWyByLFp2lriyZc9BJj2Zx+qdhdw8tg/9krYqDESOU3zuZiEJKXd1Pl994F22F5Tw2DWD+f45p5Cko45FjpuWECTmuTt/zV3Pva+v5tSMFjw8IZtu6c2iXZZIzFEgSEwrLq1gyj+WMmvZDi4ccBJ3j++vvYhEwqS/HIlZh7cXrNlZyC++3Jfrh/fUielEakGBIDHpvXW7+d7Ti3GH6dcO4ew+7aNdkkjMUyBITHF3npi3mV+/tIJe7Zoz9apserRrHu2yROKCAkFiRllFgF/+eznPLPiYMX078KdvnEGLJtqlVKSuKBAkJuwtLuO7f1/Ego17uWHkyfzk3FNJStL2ApG6pECQBm/tzkK+/XgeOw6UcN/lZ3DRGZ2jXZJIXFIgSIOWuzqf7z+9hNRGyTw7aRgDu7WJdkkicUuBIA3W43M38auZH5HZsSWPXJ3NSa2bRrskkbimQJAGp6IywJ0vreDxeZsZ0zeD+y4/Q5e4FKkH+iuTBqWotILvP72Y2at3MXF4T249vy/J2ngsUi8UCNJgbC84xHXTg0ce33Xx6Vw5tHu0SxJJKAoEaRA++qSA66YvpLi0kkevGcwIHXksUu8UCBJ1s1fnc+NTi2nVtBH/nHwmp3ZsEe2SRBKSAkGi6pkFW7jtxeVkdmzBo9cMJqNlk2iXJJKwIhYIZjYFuBBoCdwPrAfuAiqAfGCCux+q8vwmwFSgK3AIuNrd8yNVn0SXu/OH19fwwOx1jDy1PQ9cMYg07UkkElURuWKamQ0GhgMjQv/2BP4GjHf3swmGw/XVhl0LbHH3kcCTwC8iUZtEX3llgJ889yEPzF7H5YO78shV2QoDkQbgqIFgZmea2TNmtsnMNpvZRjN71sy+dIzXvQBYBbxY5Tba3bcffmmCSwpVnQ3MCt2fRTBIJM4UlVZw3fSF/HPxVn48tg+/vaQ/Kcm6kqtIQ2DuXvMMsweBA8BTwDIPPdHMMoGrgFbu/r0jjH0I6AOMA7oDM4FMd3czGw/8DDjb3Q9WGfM6cJO7f2RmKcB6d//cfodmNgmYBJCRkZGVk5MTVuNFRUWkpaWFNTZWRbvn/aUB/rSolI8LA1zbrzHDu0T+TKXR7jkaEq3nROsXat/zqFGjFrl79udmuHuNNyDjSPOONR+4D5hSZXoR0B64EXgP6FDDmGeBYaH77YDFR3t/dycrK8vDNXv27LDHxqpo9rxhV5Gf9bs3PfO2V/ytVTvr7X31c45/idave+17BvK8hs/Uoy2r7zazn5tZKoCZjTOz/zGzpFCQ7DzK2HnAGDNLMrNOQBrBb/VjgHO85o3Fc4DzQvcvCE1LHFj68X7GPziX4tJKnpk0jFGndoh2SSJSg6MFwt0E9/g5bBHQCbjnOF53BvAhMBd4AbgduAPoCLxqZrlmNtnM2pvZc6Ex04C+ZpYLfBP4zQn0IQ3UO2t28c2p79M8NZl/Tj6TM7q2jnZJInIER9u14yx3/+LhCXffZWY3EPz2f1TuHgCmVHt4xhGefllozCHg8mO9tsSOf32wjZtnLOWUjBY8fu1gOugYA5EG7WiBcKj6A+7uZhaIYD0SJ6a/t5E7Zq5gSM+2PHJ1Ni11qUuRBu9oq4z2m9mQqg+Eji8ojWxJEsvcnT/+Zw13zFzB2NMyeOK6IQoDkRhxtCWEW4CXzWwBsAXoAQwjtIpHpLpAwPnVzI94fN5mLsvqomMMRGLMEQPB3dea2SCCe/x0Az4Avu1VTjchclh5ZYApzy3lxQ8+YeLwnvz8y30x03UMRGLJEQPBzFoS3E30FXc/UH8lSawpKa/kxqcX88bKfKacdyo3jDxZYSASg462PP93YFfoX5EaFZVWcO1jC3ljZT53XtSP743qrTAQiVFH24ZQTPD0E0X1VIvEmP0Hy7j6sYUs31bAn74xgIsHdol2SSJSC0cLhAnAAODxeqpFYsiuwlImTJvPhl3FPHjlIM7t1zHaJYlILR1tldENwAfuXv2spJhZipn9IHJlSUP2yf5DfOOheWzec5Bp12QrDETixNGWEJYAr5jZUmAxUACkE9z1tC/wy8iXJw3Nlj0HueKR9yk4WM4T3x7C4B5to12SiNSRo+12OsfMzgMmAqcRPFtpPsHrHNwYOmOeJJD1u4q4cup8DpVX8tTEoXyhS+tolyQidehou51+k+AlMEcBbxG8qE1b4Bvufn/9lCcNxeodhVz5yHzcnZxJw+jbqWW0SxKROna0VUavAtsJriZ6KPRYJbA20kVJw7J8WwETps2nUXIST08aRu8OLaJdkohEwNFWGe0DckM3SVBLP97PhGnzSUtN4emJw+jRrnm0SxKRCNGVzeWIFm3exzWPLqBVs0Y8M3EYXds2i3ZJIhJBCgSpUd6mvVzz2ELS0xrz9MRhdG7dNNoliUiEKRDkcxZu2ss1jy6gQ8smPDNxGB1b6cI2IolAgSCfMX/DHq6dvpCOrYJhkKGrnIkkDJ2sXv5r/oY9XPPYQjq1akKOwkAk4UQsEMxsipnNMbOlZnZ96LFkM5tlZuOOMOZjM8sN3X4bqdrk8xZs3Mu10xdyUusmPDNpmK5/LJKAIrLKKHSpzeHACCANuMXMegOPEbzYTk1jegDL3P3LkahJjmzBxr1c89gCOrUKhUELhYFIIrJInIHCzO4AmgGZhAIBKAdKgFuBHHd/tdqY8cDtwF7gEPBjd19dw2tPAiYBZGRkZOXk5IRVY1FREWlpaWGNjVU19bx2XyX35pXQtolxy+AmtG4SX2sR9XOOf4nWL9S+51GjRi1y9+zPzXD3Or8RPLJ5NpBK8JoKq/k0fKYD42oYcxbB02Icvp93rPfJysrycM2ePTvssbGqes+LNu/1025/xUf9frbvLDgUnaIiTD/n+Jdo/brXvucjfb5G6utgCTDL3UvdfQ3Bi+y0O8aYRcDzAO7+LtDJdOmtiFn68X6unhbatVTbDESEyG1UngeMMbMkM+tEcLXR7mOMuR34CYCZDQQ2hZJM6tjhcxO1ad6YpycO1d5EIgJE7jiEGUAWMDc0PflIH+5m9iZwLnAP8HczexuoAK6LUG0JbeX2A3xr2nxaNGnE0xOH0qmVjkAWkaCIBIK7B4ApR5h3TbXpc0J39wNfiUQ9ErStKMDNj8ynaaNknpk4jC5tdG4iEflUfO1SIke0cXcx9ywsISnJeOr6oXRLVxiIyGcpEBLAx3sPcuXU9wkEnKevH0qv9om1i56IHB8FQpzbUVDClY/Mp6i0gimDm3BKhi5uIyI1UyDEsd1FpVz5yPvsKSrl8euG0K1lcrRLEpEGTIEQpwoOljNh2gK27T/Eo9cMZmC3NtEuSUQaOAVCHCoureCa6QtYn1/EQxOyGdorPdoliUgM0PUQ4kxJeSUTn8jjw60F/OWKQYzo0z7aJYlIjNASQhwprwxw49NLmLt+D7+/9AuMO71jtEsSkRiiQIgTgYDz0398yBsrd/Lri/pxyaAu0S5JRGKMAiEOuDt3zPyIF5ZsY8p5p3LVF3tEuyQRiUEKhDjwx/+s4Yl5m/nO2b24YeTJ0S5HRGKUAiHGTXt3I/e/tY7LB3fl1vMz0RnDRSRcCoQY9s9FW7nzpRWcf3pH7rq4v8JARGpFgRCj3lixk5/+80O+1DudP19+BslJCgMRqR0FQgxauGkv33t6Maef1JKHJmSTmqJTUohI7SkQYszK7Qe4bvpCOrdpymPXDiEtVccWikjdUCDEkI/3HuTqRxfQvHEKT1w3hLbNG0e7JBGJI/p6GSP2FJVy1aMLKK0I8Nx3v6irnYlInYvYEoKZTTGzOWa21MyuDz2WbGazzGxcDc83M7vPzN4xs9lmdkqkaos1xaUVXDd9IZ/sP8S0q7Ppo2saiEgERCQQzGwwMBwYEfq3p5n1BnKBfkcYNg5o5+5nA78A7o1EbbGmvDLA5KcWs2xbAQ9cMYjsHm2jXZKIxClz97p/UbM7gGZAJpAG3AKUAyXArUCOu79abcxvgeXu/lRoeqO796zhtScBkwAyMjKycnJywqqxqKiItLSGfSlJd2fqsjLmflLBtf0aM6Jro1q9Xiz0XNfUc/xLtH6h9j2PGjVqkbtnV388UtsQOgF9CH7r7w7MBDLd3Y9y8FQ6sL9qbWaW5O6Bqk9y94eBhwGys7N95MiRYRWYm5tLuGPry+9eXcXcT9bz47F9+ME5tV+DFgs91zX1HP8SrV+IXM+R2oZQAsxy91J3XwMUAe2OMWYfUHXluFcPg0Ty+NxNPJi7niuGduP7o3tHuxwRSQCRCoR5wBgzSzKzTgRXG+0+xpg5wHkAZjYCWBKh2hq8V5dv546ZHzH2tAzuvOh0nZJCROpFpFYZzQCygLmh6cl+hI0VZvYmcC4wCzjfzN4OzZoYodoatLxNe/lBzgcM7Nqa/7t8oE5JISL1JiKBEFrVM+UI866pNn1OlcnvRaKeWLF+VxHXP5FH59ZNeeTqwTRtrFNSiEj90ZHKDUR+YQlXP7qAlCTj8Wt1FLKI1D8dqdwAHCyr4NvT89hTVEbOpGF0S9dRyCJS/7SEEGUVlQG+//QSPvqkgAeuGMiArq2jXZKIJCgtIUSRu/OrmSt4c1U+d17Uj3P6ZkS7JBFJYFpCiKKpczbw5PubmXR2LyZ8sUe0yxGRBKdAiJJZy7bzm1mr+HL/jtw6LjPa5YiIKBCiYfGWfdz07AcM6taaP379DJJ0rIGINAAKhHr28d6DTHw8j4yWTZh6VTZNGulYAxFpGBQI9ajgUDnXTl9IRcB57NrBpKelRrskEZH/UiDUk/LKADc8tYjNe4p5aEIWJ7dPrNP1ikjDp91O64G7c/uLy3lv3R7uvWwAw3qlR7skEZHP0RJCPZg6ZwM5Cz/mxlG9uTSrS7TLERGpkQIhwl77aAe/fWUVF3yhEz8e2yfa5YiIHJECIYKWbyvgRzkf8IUurfnDZQO0e6mINGgKhAjZeaCE6x/Po02zRky9Kku7l4pIg6eNyhFwqKyS6x/P40BJOf/47pl0aNEk2iWJiByTAqGOBQLOT55byvJPCnh4QjanndQy2iWJiBwXrTKqY/e9uZaXl23n1nGZjD1NZy8VkdgRsSUEM5sCXAi0BO4H3gCmhd5zJfA9d6+sNuZjYH1ocp67/yxS9UXCzKWfcN+baxk/qAuTzu4V7XJERE5IRJYQzGwwMBwYEfq3J/A74G53HxF63wurjekBLHP3kaFbTIXBh1v385PnlpLdvQ2/ueR0zLRHkYjElkitMroAWAW8WOU2lOBSAsAsgkFRVRZwkpm9ZWYvm9mpEaqtzuUfKGHiE3m0S0vlbxOySE3RHkUiEnvM3ev+Rc0eAvoA44DuwEwgzd07h+aPASa4+9VVxpwFdHb3Z0P3/+zu2TW89iRgEkBGRkZWTk5OWDUWFRWRllb78wmVVTp3LyhhW1GA24Y1pWuLhrtZpq56jiXqOf4lWr9Q+55HjRq1qKbPV9y9zm/AfcCUKtOLgGI+DaBLgT9WG9MUaFRletvh5x/plpWV5eGaPXt22GMPCwQC/qOcJd79lpf8lWXba/16kVYXPcca9Rz/Eq1f99r3DOR5DZ+pkfo6Ow8YY2ZJZtYJSANeBc4Ozb8AmFNtzO3ATwDMbCCwKVR4g/XwOxt4Yck2bh7bh3Gnd4x2OSIitRKpvYxmENwmMDc0PRnYBDxiZof3Mvo3gJm9CZwL3AP83czeBiqA6yJUW52YvTqfu19dxQX9O3Hj6N7RLkdEpNYiEgjuHgCm1DBrdA3PPSd0dz/wlUjUU9fW7yriB88soW/Hlvz+si9ojyIRiQsNdwtoA1VwqJyJj+fRODmJh6/KolljHewtIvFBn2YnoDLg/DBnCVv2HuSp64fSpU2zaJckIlJnFAgn4I//WU3u6l3879dOZ6iueiYicUarjI7Tyx9u5y+z1/PNIV25cmi3aJcjIlLnFAjHYeX2A/zkuaVkdW/DHRf200ZkEYlLCoRj2H+wjElP5tGyaQoPXjlIp6UQkbilbQhHURlwvv/MEnYWlJLznWF0aKkL3YhI/FIgHMW9r69mztrd3H1JfwZ1axPtckREIkqrjI7g5Q+382Dueq4Y2o3Lh2gjsojEPwVCDdbsLGTKP5YysFtrfvnV06JdjohIvVAgVFNwqJzvPLmIZo1T+Nu3dG0DEUkcCoQqAgHn5hkf8PHeg/z1ykFkaCOyiCQQBUIVD8xexxsr87ntgr4M6dk22uWIiNQrBUJI7up8/vTGGi4e2Jmrz+wR7XJEROqdAgH4eO9BfpjzAadmtOA3F/fXkcgikpASPhBKyiuZ/NQiAu48NCGLpo21EVlEElPCH5j2P/9azvJtB5h2dTbd05tHuxwRkahJ6CWEZxduYUbeVr4/ujfn9M2IdjkiIlEVsUAwsylmNsfMlprZ9WbWw8zeNLO3zexvZpZc7flmZveZ2TtmNtvMTolUbQCbCiq5/V8fMfyUdvxoTJ9IvpWISEyISCCY2WBgODAi9G9P4HfA3e4+IvS+F1YbNg5o5+5nA78A7o1EbRA8g+kDH5TSrnlj7rt8IMlJ2ogsIhKpJYQLgFXAi1VuQ4E3QvNnEQyKqs4OPY67zwW+EInCAgHnxzOWsq/E+cuVg2jbvHEk3kZEJOZEaqNyJ6APwW/93YGZQCN399D8IqD6NSjTgf1VazOzJHcPVH2SmU0CJgFkZGSQm5t7QoVVBpzU0nIu7eUUbFhK7oYTGh7TioqKTvj/K9ap5/iXaP1C5HqOVCCUALPcvRRYY2ZFQBczs1AotAb2VBuzD2hRZdqrh0HowYeBhwGys7N95MiRJ1zcOaMhNzeXcMbGMvWcGBKt50TrFyLXc6RWGc0DxphZkpl1AtKAVwmuFoLgKqU51cbMAc4DMLMRwJII1SYiIjWI1BLCDCALmBuangxsAh4xsxRgJfBvADN7EziX4PaD883s7dCYiRGqTUREahCRQAit6plSw6zRNTz3nCqT34tEPSIicmwJfWCaiIh8SoEgIiKAAkFEREIUCCIiAigQREQkxD49eDj2mNkuYHOYw9sBu+uwnFignhNDovWcaP1C7Xvu7u7tqz8Y04FQG2aW5+7Z0a6jPqnnxJBoPSdavxC5nrXKSEREAAWCiIiEJHIgPBztAqJAPSeGROs50fqFCPWcsNsQRETksxJ5CUFERKpQIIiICBDHgWBmt5rZe6HbsGrzBprZnNDtV8czJhacaM9mlmJm00OPLTCzr0an8vCE8zMOzWtmZhvMLLN+K669MH+vp4QeW2pm19d/1bUTxu91kplNCz1/npmdEZXCa+FYn0Wh3+HFVX+H6+Tzy93j7gacBrwDGMFLeOZVm/8+0Cd0/zVg4LHGNPRbmD1fDTwQeqwdsDnafUSy3yrz/gDsBTKj3Uc9/IwHE7z2SBLQErgr2n3UQ8/jgOdCj50LvBztPuq458HAQmDH4d/huvr8itclhOHAax60meD1mVsCmFkq0Nbd14Se+0ro+UccEyPC6fmfwC9Cj33ucqUNXDj9YmaDgbbAh1GoubbC6fkCYBXwYpVbLAmn5zKguZklEbws7/Io1F0bx/osSgUuJvhzPd4xxyVeAyEd2F9luij02OF5BTXMO9qYWHDCPbt7kbsXmFkL4B/A7fVRaB054X5DV+u7m5ov3hQLwvm97kTwG+VlwHeBv5uZRbzSuhNOz3MILg2tIrh75oqIV1m3jvpZ5O7vuvvWExlzvOI1EPYR/GZwWGtgzzHmHW1MLAinZ8ysM/AG8LS7PxHxKutOOP3+FHjS3WP1vDfh9FwCzHL30tA36SKCqwdjRbg/57nu3gc4A7g7xpb2w/ksqpPPr3gNhDkE1x1iZj2Bcnc/AODuh4ACM+sV+qZ0PvDu0cbEiBPu2cw6Aa8Dv3D3R6JUd7jC+RmPA64xs1yCHxRPmNnnTvDVgIXT8zxgTGhDaycgjdg6EVw4PTcnuH4dgh+KBUBpfRdeC+F8FtXJ51dErqkcbe6+3Mxmm9kcIBmYbGYTgMbuPg24EXiC4AaYN9x9MUD1MVEqPyzh9Gxm9xFcrLzNzG4LvdT5oT+0Bi3Mn/HZh8eHQuG77r6r/qsPT5g/4w+ALGBu6GUme2grZCwIs+eNwGNmdhHBz7ifunvMBMJx9HxcY8J5bx2pLCIiQPyuMhIRkROkQBAREUCBICIiIQoEEREBFAgiIhKiQBCpJTNrZ2bHdcESM/tfM+sX6ZpEwqHdTkVqycweBP7i7sc8Z46ZpQOPu/tXIl+ZyImJywPTRCLFzG4AznL3K8zscSAPGHA4DMxsLcGjg08neBR4C4LnEtrg7pe7+x4zKzezTHdfdYS3EYkKBYLICXD3v5rZWDObTvDvZxXwUZWn9AJGAdsJnjbhS+7+kZmtN7N2ofMozQdG89mzVYpEnQJB5MTdQ/BUENnAqXz2jJt7D5+J0syK3P1wWOwDmoTu7wG61FOtIsdNG5VFTkDoHPz3Ad8GHiT44V71NMOVx/EyrYmtE8xJglAgiJyY3wP/cvdHgZeBS4EBJ/ga2cBbdV2YSG1pLyORWgrtZfSwuy85jue2AZ5w95i6frUkBi0hiNTe/wA3HOdzbwZujWAtImHTEoKIiABaQhARkRAFgoiIAAoEEREJUSCIiAigQBARkZD/B6qzq6IrCL1LAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "最大温度为63.77 C\n"
     ]
    }
   ],
   "source": [
    "lambda_ = 177\n",
    "delta = 6e-3\n",
    "L = 200e-3\n",
    "I = 800 # W/m^2\n",
    "t_water = 60\n",
    "\n",
    "# 将问题看成厚度为 2*delta，高度为 L/2 的等截面直肋的换热分析，该换热为等热流密度换热\n",
    "# 对厚度为dx的截面微元块分析，可知微分公式为：-lambda_ * delta * t''(x) = I\n",
    "\n",
    "# 试试用scipy的bvp工具求解常微分方程\n",
    "# 记y = t'(x)\n",
    "# 定义t的导数函数，返回其1阶导数和2阶导数\n",
    "N = 100\n",
    "\n",
    "\n",
    "def derivative(x, t):\n",
    "    return np.vstack((t[1], - I / (lambda_ * delta) * np.ones(len(t[1]))))\n",
    "\n",
    "\n",
    "# 定义边界条件，t(0) = t_water, t'(L/2) = 0\n",
    "# 边界条件函数的t_a和t_b参数为边界条件的起始和终止点\n",
    "def bc(t_a, t_b):\n",
    "    return np.array([t_a[0] - t_water, t_b[1]])\n",
    "\n",
    "\n",
    "x = np.linspace(0, L/2, N)\n",
    "t = np.zeros((2, x.size))\n",
    "t[0, 0] = t_water\n",
    "result = solve_bvp(derivative, bc, x, t)\n",
    "\n",
    "x_plot = np.linspace(0, L/2, N)\n",
    "t_plot = result.sol(x_plot)[0]\n",
    "plt.plot(x_plot, t_plot)\n",
    "plt.grid()\n",
    "plt.xlabel('x(m)')\n",
    "plt.ylabel('t(°C)')\n",
    "plt.show()\n",
    "print(f'最大温度为{t_plot[-1]:.2f} C')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-59"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mH = 0.2694301256218253\n",
      "Phi = 248.46 W\n"
     ]
    }
   ],
   "source": [
    "delta = 15e-3\n",
    "d_i = 120e-3\n",
    "d_o = 140e-3\n",
    "d_f = 250e-3\n",
    "lambda_ = 45\n",
    "t_i = 300\n",
    "t_f = 20\n",
    "h = 10\n",
    "\n",
    "H_p = (d_f - d_o)/2 + delta\n",
    "A_L = 2 * H_p * delta\n",
    "efficiency = fin_tip_efficiency2(H_p, A_L, lambda_, h)\n",
    "r_1 = d_o / 2\n",
    "r_2 = d_o / 2 + H_p\n",
    "area = 2 * np.pi * (r_2**2 - r_1**2)\n",
    "\n",
    "R_tube = cylindrical_wall_R(d_i/2, d_o/2, lambda_, 2*delta)\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    # 传热热阻可以看成管道导热热阻和法兰对热传热热阻的串联，流经管道和法兰的热流量相等\n",
    "    Phi, t_o = p\n",
    "    xpr1 = Phi - efficiency * h * area * (t_o - t_f)\n",
    "    xpr2 = Phi - (t_i - t_o) / R_tube\n",
    "    return xpr1, xpr2\n",
    "\n",
    "\n",
    "guess_values = [1, t_i - 10]\n",
    "Phi, t_o = root(expressions, guess_values).x\n",
    "print(f'Phi = {Phi:.2f} W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-60"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "每片肋片的热阻为：12.44 W/K\n"
     ]
    }
   ],
   "source": [
    "L = 8e-3\n",
    "# t_0 = t_L\n",
    "h = 100\n",
    "lambda_ = 200\n",
    "delta = 1e-3\n",
    "a = 100e-3\n",
    "b = 200e-3\n",
    "c = 14e-3\n",
    "\n",
    "H = L / 2\n",
    "perimeter = 2 *(a + delta)\n",
    "A_c = a * delta\n",
    "R = fin_tip_R(H, perimeter, A_c, lambda_, h)\n",
    "print(f'每片肋片的热阻为：{R:.2f} W/K')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-63"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "每米长烟道上的散热量为：673.10 W\n"
     ]
    }
   ],
   "source": [
    "t_1 = 80\n",
    "t_2 = 25\n",
    "lambda_ = 1.5\n",
    "length = 1\n",
    "\n",
    "l = 1\n",
    "d = 0.5\n",
    "\n",
    "S = 2 * np.pi * length / (np.log(1.08 * l / d))\n",
    "Phi = S * lambda_ * (t_1 - t_2)\n",
    "print(f'每米长烟道上的散热量为：{Phi:.2f} W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-65"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "热损失为：12.60 W\n"
     ]
    }
   ],
   "source": [
    "delta = 300e-3\n",
    "lambda_ = 0.8\n",
    "t_i= 400\n",
    "t_o= 50\n",
    "\n",
    "S = 0.15 * delta\n",
    "Phi = lambda_ * S * (t_i - t_o)\n",
    "print(f'热损失为：{Phi:.2f} W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-66"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "每米管道的散热量为： 7.91 W\n",
      "结冰后，每米管道的散热量为： 9.06 W\n"
     ]
    }
   ],
   "source": [
    "H_1 = 3\n",
    "H_2 = 2\n",
    "d = 500e-3\n",
    "lambda_ = 1\n",
    "t_s = 0\n",
    "t_w = 4\n",
    "length = 1\n",
    "\n",
    "if H_1 > 2 * d:\n",
    "    S_1 = 2 * np.pi * length / np.log(4 * H_1/d)\n",
    "else:\n",
    "    print('课本无对应公式')\n",
    "\n",
    "Phi_1 = lambda_ * S_1 * (t_w - t_s)\n",
    "print(f'每米管道的散热量为： {Phi_1:.2f} W')\n",
    "\n",
    "# 结冰后，地下深度为1 m的土壤层的温度跟原来相同，为0摄氏度\n",
    "if H_2 > 2 * d:\n",
    "    S_2 = 2 * np.pi * length / np.log(4 * H_2/d)\n",
    "else:\n",
    "    print('课本无对应公式')\n",
    "\n",
    "Phi_2 = lambda_ * S_2 * (t_w - t_s)\n",
    "print(f'结冰后，每米管道的散热量为： {Phi_2:.2f} W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-68"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "delta = 0.031 m\n"
     ]
    }
   ],
   "source": [
    "l_1, l_2, l_3 = 0.5, 0.75, 0.75\n",
    "lambda_ = 0.02\n",
    "t_i = -10\n",
    "t_o = 30\n",
    "Phi = 45\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    delta = p\n",
    "    S_corner = 4 * 0.15 * delta\n",
    "\n",
    "    length_1 = l_1 - 2*delta\n",
    "    length_2 = l_2 - 2*delta\n",
    "    length_3 = l_3 - 2*delta\n",
    "    length_total = 4 * length_1 + 2 * (length_2 + length_3)\n",
    "    S_edge = 0.54 * length_total\n",
    "\n",
    "    A_wall_1 = (l_1 - 2 * delta) * (l_2 - 2 * delta)\n",
    "    A_wall_2 = (l_2 - 2 * delta) * (l_3 - 2 * delta)\n",
    "    A_wall_3 = (l_3 - 2 * delta) * (l_1 - 2 * delta)\n",
    "    A_total = 2 * (A_wall_1 + A_wall_3) + A_wall_2\n",
    "    S_wall = A_total / delta\n",
    "\n",
    "    S_total = S_corner + S_edge + S_wall\n",
    "    xpr = Phi - lambda_ * S_total * (t_o - t_i)\n",
    "    return xpr\n",
    "\n",
    "\n",
    "guess_value = np.array([l_1 / 10])\n",
    "delta = root(expressions, guess_value).x[0]\n",
    "print(f'delta = {delta:.3f} m')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-73"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "芯片的工作温度为：75.34 C\n"
     ]
    }
   ],
   "source": [
    "a = 10e-3\n",
    "b = 10e-3\n",
    "delta_3 = 0.02e-3\n",
    "R_12c = 0.9e-4   # m^2-K/W\n",
    "t_oo = 25\n",
    "h = 150\n",
    "q = 1.5e4   # W/m^2\n",
    "lambda_3 = 236\n",
    "\n",
    "R_1 = 1 / (h * a * b)\n",
    "R_12 = R_12c / (a * b)\n",
    "R_3 = delta_3 / (lambda_3 * a * b)\n",
    "R_34 = 1 / (h * a * b)\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    t = p[0]\n",
    "    Phi_1 = (t - t_oo) / R_1\n",
    "    Phi_2 = (t - t_oo) / (R_12 + R_3 + R_34)\n",
    "    xpr = q * a * b - Phi_1 - Phi_2\n",
    "    return xpr\n",
    "\n",
    "\n",
    "guess_values = np.array([t_oo + 10])\n",
    "t = root(expressions, guess_values).x[0]\n",
    "print(f'芯片的工作温度为：{t:.2f} C')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-74"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Phi = 5.70 W\n"
     ]
    }
   ],
   "source": [
    "t_i = 20\n",
    "t_o = 5\n",
    "h_i = 7\n",
    "h_o = 28\n",
    "delta_1 = 12e-3\n",
    "lambda_1 = 0.16\n",
    "delta_2 = 25e-3\n",
    "lambda_2u = 0.141\n",
    "l_2u = 75e-3\n",
    "l_2d = 330e-3\n",
    "fluid_2d = 'Air'\n",
    "delta_3 = 200e-3\n",
    "lambda_3 = 0.72\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    t_12, t_23, Phi = p\n",
    "    T_12 = sc.convert_temperature(t_12, 'C', 'K')\n",
    "    T_23 = sc.convert_temperature(t_23, 'C', 'K')\n",
    "    lambda_2d = psi('L', 'T', (T_12 + T_23)/2, 'P', sc.atm, fluid_2d)\n",
    "\n",
    "    height = l_2u + l_2d\n",
    "    area = height # 单位长度面积\n",
    "    R_i = 1 / (h_i * area)\n",
    "    R_1 = delta_1 / (lambda_1 * area)\n",
    "    R_2u = delta_2 / (lambda_2u * l_2u)\n",
    "    R_2d = delta_2 / (lambda_2d * l_2d)\n",
    "    R_2 = 1 / (1/R_2u + 1/R_2d)\n",
    "    R_3 = delta_3 / (lambda_3 * area)\n",
    "    R_o = 1 / (h_o * area)\n",
    "\n",
    "    R_total = R_i + R_1 + R_2 + R_3 + R_o\n",
    "    # print(R_i*area, R_1*area, R_2*area, R_3*area, R_o*area, R_total*area)\n",
    "    xpr1 = Phi - (t_i - t_o) / R_total\n",
    "    xpr2 = t_i - t_12 - Phi * (R_i + R_1)\n",
    "    xpr3 = t_i - t_23 - Phi * (R_i + R_1 + R_2)\n",
    "    return [xpr1, xpr2, xpr3]\n",
    "\n",
    "\n",
    "guess_values = [t_i - 1, t_i - 2, 1]\n",
    "t_1, t_2, Phi = root(expressions, guess_values).x\n",
    "print(f'Phi = {Phi:.2f} W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-75"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "q_r = 37298.92 W/m\n"
     ]
    }
   ],
   "source": [
    "t_s2 = 25\n",
    "t_s1 = 150\n",
    "lambda_ = 15\n",
    "r_2 = 35e-3\n",
    "r_3 = 48e-3\n",
    "\n",
    "R = cylindrical_wall_R(r_2, r_3, lambda_, 1)\n",
    "q_r = (t_s1 - t_s2) / R\n",
    "print(f'q_r = {q_r:.2f} W/m')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-76"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "苹果中心温度为5.11 C\n",
      "苹果表面温度为5.09 C\n"
     ]
    }
   ],
   "source": [
    "Q = 4000  # J/kg\n",
    "time = 1 * sc.day\n",
    "rho = 840\n",
    "lambda_ = 0.5\n",
    "h = 6\n",
    "d = 80e-3\n",
    "t_oo = 5\n",
    "v = 0.6\n",
    "\n",
    "Q_v = Q * rho\n",
    "Phi = Q_v / time\n",
    "\n",
    "Delta_T = sphere_Delta_T_with_heat_source(d/2, lambda_, Phi)\n",
    "surface_area = np.pi * d**2\n",
    "V = 4 / 3 * np.pi * (d/2)**3\n",
    "t_r = t_oo + Phi * V / (h * surface_area)\n",
    "t_0 = t_r + Delta_T\n",
    "print(f'苹果中心温度为{t_0:.2f} C')\n",
    "print(f'苹果表面温度为{t_r:.2f} C')\n",
    "\n",
    "# # 也可以用scipy的bvp工具求解\n",
    "# # 记y = t'(r)\n",
    "# # 定义t的导数函数，返回其1阶导数和2阶导数\n",
    "# N = 100\n",
    "#\n",
    "#\n",
    "# def derivative(r, t):\n",
    "#     return np.vstack((t[1], - 2/r * t[1] - Phi/lambda_))\n",
    "#\n",
    "#\n",
    "# # 定义边界条件，t(0) = t_water, t'(L/2) = 0\n",
    "# # 边界条件函数的t_a和t_b参数为边界条件的起始和终止点\n",
    "# def bc(t_a, t_b):\n",
    "#     return np.array([t_a[1], lambda_ * t_b[1] + h * (t_b[0] - t_oo)])\n",
    "#\n",
    "#\n",
    "# r = np.linspace(0.0001, d/2, N)\n",
    "# t = np.zeros((2, r.size))\n",
    "# t[0, 0] = t_oo\n",
    "# result = solve_bvp(derivative, bc, r, t)\n",
    "#\n",
    "# x_plot = np.linspace(0.0001, d/2, N)\n",
    "# t_plot = result.sol(x_plot)[0]\n",
    "# plt.plot(x_plot, t_plot)\n",
    "# plt.grid()\n",
    "# plt.xlabel('r(m)')\n",
    "# plt.ylabel('t(°C)')\n",
    "# plt.show()\n",
    "# print(f'苹果中心温度为{t_plot[0]:.2f} C')\n",
    "# print(f'苹果表面温度为{t_plot[-1]:.2f} C')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-78"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "由于肌肉运动所造成的最大温升为1.35 C\n"
     ]
    }
   ],
   "source": [
    "r = 2e-2\n",
    "Phi_dot = 5650\n",
    "t_o = 37\n",
    "lambda_ = 0.42\n",
    "\n",
    "Delta_T = cylinder_Delta_T_with_heat_source(r, lambda_, Phi_dot)\n",
    "print(f'由于肌肉运动所造成的最大温升为{Delta_T:.2f} C')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-79"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "手柄传递的热流量为: 2.39 W\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEECAYAAAArlo9mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjtElEQVR4nO3deXxU9b3/8ddnZrKvJIGwhC3sIIsQBVEgcbe0atXqtV6s9vZy1bb6s621t+1te297b7He1mtbq9JarbUV27pUq6KihMUFZVUUEARB1rAGkpB1vr8/ZjCIIUDIycnMvJ+Pxzxm5pw5+X7yfUzec/I953zHnHOIiEj8C/hdgIiIdAwFvohIglDgi4gkCAW+iEiCUOCLiCSIkN8FHE1BQYHr169fm7evrq4mIyOj/QqKUeqHZuqLZuqLZvHWF0uWLNnlnOva0rpOG/j9+vVj8eLFbd6+vLyc0tLS9isoRqkfmqkvmqkvmsVbX5jZxqOt05COiEiCUOCLiCQIBb6ISIJQ4IuIJAgFvohIglDgi4gkCE9OyzSzfwGmRZ8mAyVAb+C3QBawE7jeOVftRfsiIvJpnuzhO+cecM6VOudKgTeAbwHfAJ5xzpUBy4HpXrS9r6aeu+esZeP+Ji9+vIhIzDIv58M3s1OB+4HxwGvAFc65LWY2BviBc+6yI14/negHQWFh4bhZs2adcJs1DY6vvVLDeUWOq0dknuyvEPOqqqrIzFQ/gPricOqLZvHWF2VlZUuccyUtrfP6StvvAT9yzjkzywf2RZdXAflHvtg5NxOYCVBSUuLaevXb2HWv8f7eyri6eq6t4u0qwpOhvmimvmiWSH3h2UFbM8sGTgdeiC7aS2T8HiAX2O1V25MGdeXD/WH2VNd71YSISMzx8iyds4G5zrlDg+kLgAuij6dGn3ti8uACHLBw3S6vmhARiTleBn4xsOqw5zOAK82sHBhBdOjGC6OKcslIggXv7/SqCRGRmOPZGL5z7hdHPN9FZM/ec8GAMTw/yPy1O3HOYWYd0ayISKcWtxdenZIfZMf+OtZWVPldiohIpxC/gV8QBGC+hnVERIA4Dvz8tAADumYwf60O3IqIQBwHPsDkwV1ZtH43tQ266lZEJL4Df1BX6hrDvPXhHr9LERHxXVwH/vjiPJKDAcrXaBxfRCSuAz89OcT44jzmrqnwuxQREd/FdeADnD20G+t3VrNxt2ZiFpHEFveBXzakGwBzV2svX0QSW9wHfr+CDIoLMnhF4/gikuDiPvAByoZ24431u6mpb/S7FBER3yRE4J89tBv1jWFeXefZjMwiIp1eQgT+af3yyEgO6mwdEUloCRH4yaEAZw0qYO7qCrz8SkcRkc4sIQIfIsM62yprWb39gN+liIj4ImECvzR6euYrOj1TRBJUwgR+YXYqI3vlMGfVDr9LERHxRcIEPsB5wwtZ/tE+Kg7U+l2KiEiHS7jAdw5eXqVhHRFJPAkV+EO7Z1HUJY2X3tOwjogknoQKfDPjvOGFLFy3i+o6XXUrIokloQIfIsM69Y1hFqzV3DoiklgSLvBP75dHTloSL76rYR0RSSwJF/ihYIBzhnbjlTUVNDaF/S5HRKTDJFzgQ2RYZ19NA299uNfvUkREOkxCBv7kwV1JDgV48b3tfpciItJhEjLwM1JCTB5UwAsrt2syNRFJGAkZ+AAXntKDrZW1rNhc6XcpIiIdwrPAN7PbzGyBma0ws6+YWT8ze9nM5pnZfWYW9Krt43HesEJCAeP5ldv8LENEpMN4EvhmdhowCZgSve8P3AHMcM5NibZ7sRdtH6+c9CTOHFjA8+9oWEdEEoN5EXZm9iMgHRgKZAK3A38F+jvnnJldCkx2zn3jiO2mA9MBCgsLx82aNavNNVRVVZGZmdnqa+ZtbuDBlfX858RU+mb7+g+HZ46nHxKF+qKZ+qJZvPVFWVnZEudcSUvrQh612QMYDFwI9AWeAZJc86dLFZB/5EbOuZnATICSkhJXWlra5gLKy8s51vajqut5+L05VKT04kulQ9vcVmd2PP2QKNQXzdQXzRKpL7waw68FnnPO1Tnn3icS8LlmZtH1uYDv3yiel5HMGcX5PKdhHRFJAF4F/uvAuWYWMLMeRIZ1ZgOTo+unAgs8avuEXHhKdzbsqmbNDn31oYjEN68C/y/A28BrwJPAjcBtwA/NbD5QDzztUdsn5IIR3TGD597RRVgiEt88GcN3zoWJBPyRzvaivZPRNSuF8f3zePbtrdx67iCaR51EROJLwl54dbjPje7JBzureW/bfr9LERHxjAIfuOiUHoQCxtMrtvpdioiIZxT4RM7WOWtQAf9YsY1wWGfriEh8UuBHXTy6J1v2HWTpJk2ZLCLxSYEfdf6I7qSEAjyjYR0RiVMK/KjMlBDnDOvGs+9s0zdhiUhcUuAf5uLRPdlVVc/r632/CFhEpN0p8A9TOqQbmSkhnl6uYR0RiT8K/MOkJgW5YER3nl+5ndqGJr/LERFpVwr8I1w+thdVdY28+N4Ov0sREWlXCvwjTCjOp2dOKk8s3ex3KSIi7UqBf4RAwPj82F7Mf38nFQdq/S5HRKTdKPBb8PlTiwg7dPBWROKKAr8FA7tlMrp3Lo8v3eJ3KSIi7UaBfxSXj+3Fqm37eW+rZtAUkfigwD+Kz43qSVLQdPBWROKGAv8oumQkc/bQbjy1fAsNmmpBROKAAr8VV5b0ZldVPa+srvC7FBGRk6bAb8WUwV3plpXCY2995HcpIiInTYHfilAwwBdKiihfU8H2Sp2TLyKxTYF/DFeW9Cbs4G9LtJcvIrFNgX8MffMzOKM4n8cWf6SvPxSRmKbAPw7/dHpvPtpzUPPki0hMU+AfhwtGdCcnLYlZOngrIjFMgX8cUpOCfP7UXrywcju7q+r8LkdEpE0U+Mfpi+P7UN8U5q9LdOWtiMQmBf5xGlyYxfj+efxp0UYdvBWRmKTAPwHTzujLR3sOMm/tTr9LERE5YZ4Fvpm9ZWbl0duDZtbPzF42s3lmdp+ZBb1q2yvnD+9OQWYKj7y+0e9SREROmCeBb2bJQMg5Vxq9XQ/cAcxwzk2JtnuxF217KTkU4OrTe/PKmgo+2lPjdzkiIifEqz38kUCGmb0U3aufAIwH5kTXPwdM8qhtT119eh8MePTNTX6XIiJyQsy59j8AaWYjgLOAmcAgYDaQ4pzrFV1/LjDNOfelI7abDkwHKCwsHDdr1qw211BVVUVmZmabt2/N3UtrWbeviV+UppMUME/aaC9e9kOsUV80U180i7e+KCsrW+KcK2lpXcijNtcC61zk0+R9M9sFjDUziy7LBT512apzbiaRDwlKSkpcaWlpmwsoLy/nZLZvTajXLv75gUVU5gziinFFnrTRXrzsh1ijvmimvmiWSH3h1ZDOdcAvAcysF5ANPA1Mjq6fCizwqG3PnTkwn8GFmfx+4Qa8+A9JRMQLXgX+Q0C6mS0EZgHXA98Cfmhm84F6Ih8AMcnM+PKZ/Xlv234WbdjjdzkiIsfFkyEd51w9MK2FVWd70Z4fLj21F3fMXs3vF25gQnG+3+WIiByTLrxqo9SkIF8c34eXVu1g026doikinZ8C/yRMm9CPoBkPvfah36WIiByTAv8kdM9JZeqoHvxl8Ufsr23wuxwRkVYp8E/SV84qpqqukUcX6UIsEencFPgnaWRRDhMH5PPAwg3UNTb5XY6IyFEp8NvBDVMGUHGgjr8v2+p3KSIiR6XAbweTBhUwvEc2983/QHPli0inpcBvB2bGv00pZv3Oauas2uF3OSIiLVLgt5OpI3tQ1CWN++Z9oOkWRKRTUuC3k1AwwL9OKmbppn2abkFEOiUFfju66rTeFGSm8KtX1vpdiojIpyjw21FqUpB/m1zMq+t2s2Sj9vJFpHNR4Lezayb0IS8jmV++vM7vUkREPkGB387Sk0N8ZVJ/5r2/k+Uf7fO7HBGRj7Ua+GY20cweNbMPzWyjmW0ws8fM7MyOKjAWXXtGP3LTk/i1xvJFpBM5auCb2b3AJcBPgf7Oub7Ouf7AD4GpZnZPB9UYczJTQnz5zP7MWVXByi2VfpcjIgK0vof/I+fc7c65t91hJ5Y751Y7574L/Jf35cWu687sR05aEj9/cY3fpYiIAK0H/i4z+66ZpQCY2YVm9gMzCwA453RJaSuyU5O4YcoA5q7ZyeIPdcaOiPivtcCfAfQ+7PkSoAfwM08riiNfmtiXgswU7nxhja6+FRHftRb4ZznnbnTO1QE453YCNwFndUhlcSA9OcTXygawaMMeFq7b5Xc5IpLgWgv8g0cuiI7lh70rJ/5cPb4PvXLT+F/t5YuIz1oL/H1mdvrhC8zsNKDO25LiS0ooyM3nDGTF5kpeeFeHPUTEP6FW1t0OPGtmbwKbgH7ABOALHVBXXLl8bBEz56/nZ7NXc86wbiQFdb2biHS8oyaPc24tMBZ4BtgLPAWMcM4t6ZjS4kcoGODfLxrG+l3VzHpT330rIv446h6+mWUD5wLPO+f2d1xJ8emcYd0Y3z+P/5uzlktP7UVWapLfJYlIgmltbOERYGf0Xk6SmfG9qcPYXV3P/fPW+12OiCSg1gK/GhgMVHVQLXFvVFEuF4/uye8Wrmdb5adOghIR8VRrgT8NWA5c2zGlJIbbLhhC2MGdszXlgoh0rNYC/yZguXOu8cgVZhYys5tb+8Fmlm5m681sqJnlmdnTZjbXzP5iZhknW3is6p2Xzr9O6s8Ty7awZONev8sRkQTSWuAvA543szvN7Goz+4yZTYvOkvlidH1rfgzkRh/fDjzjnCsj8l/D9JMrO7bdVDqQwuwU/vOZdwmHdTGWiHSM1k7LXABcAKwFhgMXA4OA1cA50fUtil6glQe8HV00GXgu+vg5YNJJVx7DMlJC/PtFw3h7cyV/W7LZ73JEJEHY0S73N7OriYR8GfAKYEAQGO2cG3LUH2gWAl4ArgL+BtwAPA2c6pyrNrOBwAPOuSktbDud6N5/YWHhuFmzZrX5F6uqqiIzM7PN23vNOcd/L6qloibMjEnppCeZJ+109n7oSOqLZuqLZvHWF2VlZUuccyUtrWvtStvZwDYgH7g/uqyJyB5/a74N/NE5t8vs4xDbC2QROfMnF9jd0obOuZnATICSkhJXWlp6jKaOrry8nJPZviMUDKrk4nsW8lZtN3543ghP2oiFfugo6otm6otmidQXRw1859xeoDx6OxEXAmEzuw4YAzwMvENkeOgPwFTgqMNBiWRkUQ7XjO/DH177kMvHFnFKrxy/SxKRONbuk7o45yY750qdc6U0n9Z5O3ClmZUDI4juxQvcdsFQ8jJS+N6T79CkA7gi4iFPZ/GKBv9q59wu59zU6PMrnXPVXrYbS3LSkviPzw5jxeZK/rxoo9/liEgc07SNncDFo3ty1sACfjZ7DRUHav0uR0TilAK/EzAzfnzpKdQ1hfnR0+/6XY6IxCkFfifRvyCDW84ZxHPvbGf2ym1+lyMicUiB34lMn1zMiJ7ZfP+pd9lXU+93OSISZxT4nUhSMMDPrhjFvpp6/usf7/ldjojEGQV+JzOiZw43TBnAE0u3MHdNhd/liEgcUeB3Ql8/ZyCDumVy+9/e1tCOiLQbBX4nlBIKctdVY9hTXc/3nlrJ0eY7EhE5EQr8TuqUXjn8v3MH8ezb23h6xVa/yxGROKDA78RumDKAsX1y+Y+nVrJ1n74SUUROjgK/EwsFA/ziyjE0hh3f/MsKzbUjIidFgd/J9SvI4EcXj+D19bv5zdx1fpcjIjFMgR8DvjCuiEvG9OSuOe/z5oY9fpcjIjFKgR8DzIyfXHoKvfPSuWXWMp2qKSJtosCPEVmpSfzq6lPZVVXHt/66Ql9+LiInTIEfQ0YV5fLdzwxjzqoK7p33gd/liEiMUeDHmOsm9uNzo3vy8xfX8Oq6XX6XIyIxRIEfY8yMGZeNZEDXTL7+6DKdny8ix02BH4MyUkLcN20c9Y1hbnxkCbUNTX6XJCIxQIEfowZ0zeTnV45mxeZKvvP425pvR0SOSYEfwy4Y0Z1vnjeYp5Zv5f756/0uR0Q6uZDfBcjJ+drZA1mz4wB3zF7NoG6ZnDOs0O+SRKST0h5+jDMz7rxiNKf0zOHmR5fx7tZKv0sSkU5KgR8H0pKD/O5LJWSnJfHlh97SmTsi0iIFfpwozE7lwetPo6auiesffIv9tQ1+lyQinYwCP44M7Z7NfdPG8cHOKm744xLqGnW6pog0U+DHmTMHFnDH5aN47YPd3PrYcsI6XVNEonSWThy6fFwRe2vq+cmzq6jaG6Ks1GFmfpclIj7zZA/fzAJm9lszW2hmC8xspJn1M7OXzWyemd1nZkEv2paIr0wq5qtlA5i/uZE7Zq/xuxwR6QS8GtL5HBB0zp0F/AD4H+AOYIZzbkq03Ys9aluivnX+EMp6h7hv3gfcPWet3+WIiM88CXzn3N+B6dGn/YA3gfHAnOiy54BJXrQtzcyMacOTuXxsEXfNeZ979BWJIgnNszF851yjmT0IXAb8O5Dkmid8qQLyj9zGzKYT/aAoLCykvLy8ze1XVVWd1Pbxoqa6mqldYWuPIHe+sIZNH27gov5JfpflC70nmqkvmiVSX5jXk26ZWSGwFEgH8pxzzsyuACY6575xtO1KSkrc4sWL29xueXk5paWlbd4+Xhzqh8amMP/vseX84+1t3H7hUG4sHeB3aR1O74lm6otm8dYXZrbEOVfS0jpP9vDN7Fqgp3NuBlADHADeACYD84CpwNNetC0tCwUD/N9VYwiYccfs1RxsaOLWcwfp7B2RBOLVkM7jwMNmNh8w4BZgLfA7MwsBq1Dgd7hQMMBdV40hJRTgly+vpa6hie9cNFShL5IgPAl851w1cHkLq872oj05fsGAccflo0hJCnD//PXsr23gx5ecQiioa/BE4p0uvEpAgYDx40tOITctmV/PXcfuqnp+efWppCbp0giReKbdugRlZnzrgiH86HPDeWnVDq79/ZtU1mjCNZF4psBPcNed2Z+7/+lUlm3ay2X3vsqm3TV+lyQiHlHgCxeP7skf/2U8u6rq+fxvXmXppr1+lyQiHlDgCwATivN54qaJZKaGuHrmG/x9+Ra/SxKRdqbAl48N6JrJkzedyeiiXG6ZtZyfPreKprCmVxaJFwp8+YS8jGQe+cp4/nlCH+6fv57rH3pLB3NF4oQCXz4lORTgJ5eO5H8+P5LXP9jF1F8t4J3N+nJ0kVinwJej+uL4Pjz2b2cQDjsuv/c1HnljI17PvSQi3lHgS6vG9unCP26exIQB+Xz/qZV8/dFlVB7UEI9ILFLgyzHlZSTz0HWncdsFQ3h+5XY+c/cClmzc43dZInKCFPhyXAIB46tlA/nrDWcQCMCV97/BXS+9T0NT2O/SROQ4KfDlhIzt04Vnb57ExaN7cvfLa7nsN6+xdscBv8sSkeOgwJcTlp2axF1XjeHea8ayZd9Bpv5qIffN+4BG7e2LdGoKfGmzi0b24MVbJ1M2pCsznl/NJfe8ysotOn1TpLNS4MtJKchM4b5/Hse914yl4kAdl9zzKj99bhXVdY1+lyYiR1Dgy0kzMy4a2YM5t07hC+OKuH/+es79xTyef2ebztsX6UQU+NJuctKTmHH5KB6/8Qxy0pK48U9Lufb3b+qgrkgnocCXdjeubx7/+PpZ/OCzw1n+0T4uvHsBP/z7SvZW1/tdmkhCU+CLJ0LBAF8+qz/zbivj6tN788c3NjLlzrncN+8Dahua/C5PJCEp8MVTeRnJ/OTSkTx/y2TG9e3CjOdXU/a/5fzlrY90GqdIB1PgS4cY0j2LB68/nUf/dQLdslP59uNvc+4v5vHkss2ac1+kgyjwpUOdMSCfp26ayMxp40hNCnLrYys4/655PLF0s/b4RTymwJcOZ2acP6I7z908id9cM5akYIBv/GUFZ/98Hn9etElj/CIeUeCLbwIB4zMje/DczZP47bUldMlI5rtPvsNZd7zCr15eq7N6RNpZyO8CRAIB47zhhZw7rBuvr9/NzPnr+flL73NP+TouG1vEdRP7Mbgwy+8yRWKeAl86DTNj4oACJg4oYM32A/x+4QYeX7KZPy/axMQB+Uyb0JdzhxeSFNQ/piJtocCXTmlI9yzuuGIU37loKLPe+ohH3tjIjX9aStesFK4q6c2VJb3pk5/ud5kiMUWBL51al4xkbiwdwPTJxcx7v4I/vbGJe8rX8eu56zijOJ8rTyvighHdSU/WW1nkWDz5KzGzEPA7YACQAvwYeAd4INrmKuCrzjmdjiHHJRgwzh5ayNlDC9m67yCPL9nMX5ds5tbHVpCevJILR3Tn82N7cUZxPiEN+Yi0yKvdomuAKufcJDMrAJYAbwAznHMvmdlM4GLgSY/alzjWMzeNr58ziK+WDeTND/fw1LItPPvONp5YtoWCzGQuOqUHnx3Vg9P65REImN/linQa5sX0tWaWCQSdc5VmlgcsAxzQ3znnzOxSYLJz7htHbDcdmA5QWFg4btasWW2uoaqqiszMzDZvHy8SpR/qmxwrdjaxaFsjK3Y20RCGnBRjXLcg4wpDDMkLUFtTnRB9cTwS5X1xPOKtL8rKypY450paWudJ4H/8w82ygL8DDwE/dc71ii4/F5jmnPvS0bYtKSlxixcvbnPb5eXllJaWtnn7eJGI/VBd18icVTt44d3tzF29k4MNTWSnhhjexfHF0lFMGdyVnLQkv8v0VSK+L44m3vrCzI4a+J4d6TKzXsATwG+dcw+b2X+ZmbnIJ0wusNurtiWxZaSEuGRMLy4Z04uD9U3MX7uTOe/tYPY7m7n50WUEA8bYPrmUDunGlMFdGd4jW0M/khC8OmjbA3gRuMU5Nye6eDEwGZgHTAWe9qJtkcOlJQe5YER3LhjRnYsK9pBTPIbyNRXMXVPBnS+s4c4X1pCfkczEgQWcNTCfM4oL6J2Xhpk+ACT+eLWH/x0gH/i+mX0/uuyrwK+iZ/CsQoEvHSxgxri+XRjXtwvfPH8IFftrWbhuFwvW7mLhul08s2IrAL1y05hQnM/4/nmc1j+Pfvnp+gCQuOBJ4DvnbgFuaWHV2V60J9IW3bJTuWxsEZeNLcI5x7qKKl5fv5vX1u1m7poKHl+6GYCuWSmM6xP5oBjbN5cRPXNITQr6XL3IidPVKiJEpnUYVJjFoMIsrj2jH845PthZxaINe3hrwx6WbtrH7He3AxAKGMN6ZDOmdy4ji3IYVZTDwK6ZOv9fOj0FvkgLzIyB3bIY2C2La8b3BaDiQC3LNu1j+Uf7WL5pH08u28If39gIQGpSgKHdsxnRM5sRPXMY1iOLwYVZZKToT0w6D70bRY5Tt6zUjw8AA4TDjg27q3lncyVvb67k3a2VPL1iK39atAkAM+iTl87gwiwGF2YyuDCLgd0yGdA1U0NC4gsFvkgbBQLGgK6RAL/01F4AOOfYvPcgq7cfYPW2/azavp/3d1TxyuqKj7/K0Qx6d0mnuGsGxQWZ9O+aQXFBBn3z0+mZk6ZTRMUzCnyRdmRm9M5Lp3deOucNL/x4eV1jExt2VbOuourj2/qd1Sxav4eDh33DV3IoQO8uafTJS6dvfgZFXdLonZdOUZc0irqkJ/wFY3JyFPgiHSAlFGRo92yGds/+xPJw2LF9fy0f7q5m4+4aPtwVud+0p4a3PtxLVV3jJ16flRKiZ24aPXNT6ZGbRs+cVLrnpNE9O5XuOSkUZqeSmRLSaaTSIgW+iI8CAYsGeBoTB3xynXOOfTUNbN57kI/21rBl70G27Ivctu47yIrNlexp4Wsg05ODFGan0jUrJXLLjNwXZCZTkJlCQWYKO2vCHKxvIi1ZxxISiQJfpJMyM7pkJNMlI5mRRTktvqa2oYntlbVs31/Ljv21bK+speJAHTv2R+5Xbd3P/AN1HDjiPwWA2+bPJi0pSF5GMl0ykuiSnhy9JZGbnkxuelLklpZMdloSOWlJZKeFyE5N0kHnGKXAF4lhqUlB+hVk0K8go9XXHaxvYldVHbuq6thdVc9rS9+ma1Exu6vq2FvTwJ7qOvbUNLBpTw17q+vZX/vpD4jDJYcCZKcmkZ0aIjM1RFZqiKyUJDJTQ2SmRG7pKUEyU0JkJIfISAmSkRIiPTlIenLkPi05SFpS5HlQB6o7hAJfJAGkJQc/PpgMEKpIorR0wFFf39gU5kBtI/sONrC3pp79BxuoPNjA/oMN7K9tZH9tA/sPNnKgtoEDtZH7nQfqqK5r4kBtA1V1jYRPYCLe5FCAtKTIB0BqUoDUpGD0Fn0cCpKSFCAlFCAlFCQlFCA5+jg5+jg5FCAlGCApZCQHgyQFjaRggKRggNDHj41QIHofDBAKGHtrw+w8UEcoYASDRihgBCxyHwxYXB0PUeCLyKeEgoGPh5P60/p/Dy1xzlHbEKaqrpHqukZq6puoqW+kur6Jg/WNVNc1UdMQeVxT38TB+iZqG5qobQhT03DocRN1DWH2VNdT1xCmtjHyvL4pTG1DE/WNYRpP5FOlNeVzjroqYJFvXAtY5AMgaEYgYAQsMj+TWfPjgEWG4gIBMAwzMCLLDCD6/JAjP0wOTVc/rEc2v/7i2Pb53Q6jwBeRdmdmkSGb5CBds1I8a6cp7KhvDEduTZFbw6HHjWEamsI0NLnofZjGJkdjOLLs0P17q1YzYNBgGpvCNIUdTWFHY/T+45tzhA977ByEXeS5IxLUTeFDyyPrnDu0jo9f84mPJ3f4Q4cd+igw6Juf7kl/KfBFJGYFA80fLG1VXvUBpRP6tmNVnZdmexIRSRAKfBGRBKHAFxFJEAp8EZEEocAXEUkQCnwRkQShwBcRSRAKfBGRBGGHLuXtbMxsJ7DxJH5EAbCrncqJZeqHZuqLZuqLZvHWF32dc11bWtFpA/9kmdli51yJ33X4Tf3QTH3RTH3RLJH6QkM6IiIJQoEvIpIg4jnwZ/pdQCehfmimvmimvmiWMH0Rt2P4IiLySfG8hy8iIodR4IuIJIiYDHwz+46ZvRq9TThi3almtiB6+8/j2SaWnWhfmFnIzB6KLnvTzD7nT+Xtry3vi+i6dDNbb2ZDO7Zi77Txb+S26LIVZvaVjq+6/bXh7yNgZg9EX/+6mY3xpXCvuENfxRUjN2A4MJ/IV0P2BRYfsf4NYHD08QvAqcfaJlZvbeyLLwG/ji4rADb6/Xv41ReHrfs5sAcY6vfv4eP74jTgaSI7gdnAf/v9e/jUDxcCf40uOx941u/foz1vsbiHPwl4wUVsBEJmlg1gZilAnnPu/ehrn4++/qjbxLi29MXjwPeiy8IdXbCH2tIXmNlpQB7wtg81e6UtfTEVWA08ddgt1rWlH+qBDDMLAFnASh/q9kwsBn4+sO+w51XRZYfWVbawrrVtYtkJ94Vzrso5V2lmWcDfgP/oiEI7wAn3hZmFgBnAbR1RYAdqy99IDyJ7+V8AbgAeMTPzvFJvtaUfFhD5D2c1kdM13/O8yg4Ui4G/l8gn7yG5wO5jrGttm1jWlr7AzHoBc4A/O+ce9rzKjtGWvvg28EfnXDzNowJt64ta4DnnXF10r7eKyJBfLGvre+I159xgYAwwI05GA4DYDPwFRMbWMLP+QINzbj+Ac+4gUGlmxdG9k4uAha1tE+NOuC/MrAfwIvA959zvfKrbC215X1wIXGdm5UT+uB82sxYnnYoxbemL14FzowctewCZxP6EYm3phwxge3T73UT+C6jr6MK9EvK7gBPlnFtpZnPNbAEQBG40s2lAsnPuAeBrwMNEDtTMcc4tBThyG5/Kb1dt6Qszu5vIv67fN7PvR3/URdE/gJjVxvfF5EPbR0P/Bufczo6vvn218X2xHBgHvBb9MTe66JHLWNXGftgAPGhmlxDJx2875+Im8HWlrYhIgojFIR0REWkDBb6ISIJQ4IuIJAgFvohIglDgi4gkCAW+yDGYWYGZHdeXZJjZT8xshNc1ibSFTssUOQYzuxe4xzl3zHlVzCwf+INz7rPeVyZyYmLuwisRL5nZTcBZzrkvmtkfgMXA6ENhb2ZriVyVegqRK5aziMxBs94590/Oud1m1mBmQ51zq336NURapMAXOYxz7jdmdp6ZPUTk72M18O5hLykGyoBtRC69P9M5966ZfWBmBdF5eRYBZ0e3Fek0FPgin/YzIlMMlABD+OSsinucc5sBzKzKOXfow2AvkBp9vBso6qBaRY6bDtqKHCY6T/rdwL8A9xIJ78On0m46jh+TS+xPPCZxSIEv8kl3An93zv0eeBa4Ahh9gj+jBHilvQsTOVk6S0fkGKJn6cx0zi07jtd2AR52zsXNdwVL/NAevsix/QC46Thf+03gOx7WItJm2sMXEUkQ2sMXEUkQCnwRkQShwBcRSRAKfBGRBKHAFxFJEP8fRfgJDP/OmlEAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "温度为35°C的位置离手柄与锅体相接部分0.017 m\n"
     ]
    }
   ],
   "source": [
    "h_i, h_o = 2, 10\n",
    "d_i, d_o = 25e-3, 30e-3\n",
    "H = 90e-3\n",
    "t_f = 15\n",
    "t_w = 70\n",
    "lambda_ = 1.3  # 缺少条件，课本并未给出，这里取黏土的导热系数\n",
    "t_search = 35\n",
    "\n",
    "r_i, r_o = d_i / 2, d_o / 2\n",
    "delta = r_o - r_i\n",
    "H_p = H + delta/2\n",
    "\n",
    "perimeter = np.pi * d_i + np.pi * d_o\n",
    "A_c = np.pi * (r_o**2 - r_i**2)\n",
    "efficiency = fin_tip_efficiency(H_p, perimeter, A_c, lambda_, 1)\n",
    "area_i = np.pi * d_i * H_p\n",
    "area_o = np.pi * d_o * H_p\n",
    "Phi = efficiency * (h_i * area_i + h_o * area_o) * (t_w - t_f)\n",
    "print(f'手柄传递的热流量为: {Phi:.2f} W')\n",
    "\n",
    "# # 用scipy的bvp工具求解具体的温度分布\n",
    "# 记y = t'(x)\n",
    "# 定义t的导数函数，返回其1阶导数和2阶导数\n",
    "N = 100\n",
    "\n",
    "\n",
    "def derivative(x, t):\n",
    "    k = 4 * (d_o*h_o + d_i*h_i) / (lambda_ * (d_o**2 - d_i**2))\n",
    "    return np.vstack((t[1], k * (t[0] - t_f)))\n",
    "\n",
    "\n",
    "# 定义边界条件，t(0) = t_water, t'(L/2) = 0\n",
    "# 边界条件函数的t_a和t_b参数为边界条件的起始和终止点\n",
    "def bc(t_a, t_b):\n",
    "    return np.array([t_a[0] - t_w, t_b[1]])\n",
    "\n",
    "\n",
    "x = np.linspace(0, H, N)\n",
    "t = np.zeros((2, x.size))\n",
    "t[0, 0] = t_w\n",
    "result = solve_bvp(derivative, bc, x, t)\n",
    "\n",
    "x_plot = np.linspace(0, H_p, N)\n",
    "t_plot = result.sol(x_plot)[0]\n",
    "plt.plot(x_plot, t_plot)\n",
    "plt.grid()\n",
    "plt.xlabel('x(m)')\n",
    "plt.ylabel('t(°C)')\n",
    "plt.show()\n",
    "\n",
    "arg = find_nearest(t_plot, t_search)\n",
    "print(f'温度为{t_search}°C的位置离手柄与锅体相接部分{x_plot[arg]:.3f} m')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-80"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "小屋内的空气平均温度为：-10.73 C\n"
     ]
    }
   ],
   "source": [
    "r = 1.8\n",
    "delta = 0.5\n",
    "lambda_ = 0.15\n",
    "t_oo = -40\n",
    "h_o = 15\n",
    "h_i = 6\n",
    "t_ground = -20\n",
    "Phi = 960\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    t = p[0]\n",
    "    area_ground = np.pi * (r**2)\n",
    "    Phi_ground = h_i * area_ground * (t - t_ground)\n",
    "\n",
    "    r_m = r + delta / 2\n",
    "\n",
    "    R_2 = spherical_wall_R(r, r+delta, lambda_)\n",
    "    R_1 = 1 / (h_i * (2 * np.pi * r**2))\n",
    "    R_3 = 1 / (h_o * (2 * np.pi * (r+delta)**2))\n",
    "\n",
    "    R_roof = R_1 + R_2 + R_3\n",
    "    Phi_roof = (t - t_oo) / R_roof\n",
    "    xpr = Phi_ground + Phi_roof - Phi\n",
    "    return xpr\n",
    "\n",
    "\n",
    "guess_value = [t_ground+1]\n",
    "t = root(expressions, guess_value).x[0]\n",
    "print(f'小屋内的空气平均温度为：{t:.2f} C')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-81"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "外边面的温度为：569.14 C\n"
     ]
    }
   ],
   "source": [
    "T_av = 470\n",
    "q = 2500\n",
    "t_i = 65\n",
    "lambda_ = [0.047, 0.012, 0.038]\n",
    "delta = [0.8e-3, 0.55e-3, 3.5e-3]\n",
    "delta_12 = delta_34 = 1e-3\n",
    "h_rad = 4 * sc.sigma * T_av**3\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    t_2, t_4, R_total = p\n",
    "    R_a = delta[0] / lambda_[0]\n",
    "    R_b = delta[1] / lambda_[1]\n",
    "    R_c = delta[2] / lambda_[2]\n",
    "    t_1 = t_i + q * R_a\n",
    "    t_12 = (t_1 + t_2)/2\n",
    "    lambda_12 = psi('L', 'T', t_12, 'P', sc.atm, 'Air')\n",
    "    R_12cond = delta_12 / lambda_12\n",
    "    R_12 = 1 / (h_rad + 1 / R_12cond)\n",
    "    xpr1 = t_2 - t_1 - q * R_12\n",
    "    t_3 = t_2 + q * R_b\n",
    "    t_34 = (t_3 + t_4)/2\n",
    "    lambda_34 = psi('L', 'T', t_34, 'P', sc.atm, 'Air')\n",
    "    R_34cond = delta_34 / lambda_34\n",
    "    R_34 = 1 / (h_rad + 1 / R_34cond)\n",
    "    xpr2 = t_4 - t_3 - q * R_34\n",
    "    xpr3 = R_total - (R_a + R_12 + R_b + R_34 + R_c)\n",
    "    return xpr1, xpr2, xpr3\n",
    "\n",
    "\n",
    "guess_values = [t_i+10, t_i+20, 1]\n",
    "t_2, t_4, R_total = root(expressions, guess_values).x\n",
    "t_o = t_i + q * R_total\n",
    "print(f'外边面的温度为：{t_o:.2f} C')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-86"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mH = 0.2740640638812596\n",
      "肋片的效率为: 97.57%\n",
      "肋面的总效率为: 99.96%\n",
      "该热沉能散发的热量为：19.74 W\n"
     ]
    }
   ],
   "source": [
    "t_0 = 75\n",
    "t_oo = 25\n",
    "h = 20\n",
    "delta = 2e-3\n",
    "H = 25e-3\n",
    "W = 12e-2\n",
    "L = 18e-2\n",
    "lambda_ = 180\n",
    "number = 6\n",
    "\n",
    "H_p = H + delta/2\n",
    "eta_fin = fin_tip_efficiency2(H_p, H_p*delta, lambda_, h)\n",
    "print(f'肋片的效率为: {eta_fin:.2%}')\n",
    "A_f = number * H_p * delta\n",
    "A_r = (W - number*delta) * L\n",
    "eta_overall = overall_fin_surface_efficiency(A_f, A_r, eta_fin)\n",
    "print(f'肋面的总效率为: {eta_overall:.2%}')\n",
    "R_overall = 1 / (h * (A_f + A_r) * eta_overall)\n",
    "# 缺少条件，未知底层厚度，这里取肋片的厚度delta\n",
    "R_conduction = delta / (lambda_ * W*L)\n",
    "Phi = (t_0 - t_oo) / (R_overall + R_conduction)\n",
    "print(f'该热沉能散发的热量为：{Phi:.2f} W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-87"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "需要3个针肋\n"
     ]
    }
   ],
   "source": [
    "t_w = 70\n",
    "l_1, l_2 = 10e-2, 16e-2\n",
    "H, d = 3e-2, 4.2e-2\n",
    "lambda_ = 15\n",
    "t_f = 20\n",
    "h = 70\n",
    "Phi = 80\n",
    "\n",
    "perimeter = np.pi * d\n",
    "A_c = np.pi * d**2 / 4\n",
    "H_p = H + d/4  # H_p = H + A_c/perimeter\n",
    "R_fin = fin_tip_R(H_p, perimeter, A_c, lambda_, h)\n",
    "\n",
    "\n",
    "def expressions(p):\n",
    "    number = p\n",
    "    R_fins = R_fin / number\n",
    "    R_base = 1 / (h * (l_1*l_2 - number*np.pi*d**2/4))\n",
    "    xpr = Phi - (t_w - t_f) * (1/R_base + 1/R_fins)\n",
    "    return xpr\n",
    "\n",
    "\n",
    "guess_values = 1\n",
    "number = root(expressions, guess_values).x[0]\n",
    "print(f'需要{int(np.ceil(number))}个针肋')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-91"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "接触热阻为：1.73e-04 m^2-K/W\n"
     ]
    }
   ],
   "source": [
    "y_list = np.array([51.2, 53.3, 54.4, 56.9, 59.5, 61.6])\n",
    "x_list = np.array([5, 25, 40, 60, 75, 95])*1e-3\n",
    "\n",
    "lambda_ = 42.8\n",
    "x = 50e-3   # 题目错误，应该是50e-3\n",
    "result1 = np.polyfit(x_list[:3], y_list[:3], 1)\n",
    "t_50_l = np.poly1d(result1)(x)\n",
    "result2 = np.polyfit(x_list[3:], y_list[3:], 1)\n",
    "t_50_r = np.poly1d(result2)(x)\n",
    "Delta_T = t_50_r - t_50_l\n",
    "\n",
    "q = lambda_ * (result1[1] + result2[1]) / 2\n",
    "R_c = Delta_T / q\n",
    "print(f'接触热阻为：{R_c:.2e} m^2-K/W')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-92"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "当空气流速为0.6 m/s时，散热量为：272.92 W\n",
      "当空气流速为1.0 m/s时，散热量为：370.73 W\n",
      "当空气流速为1.5 m/s时，散热量为：469.95 W\n",
      "当空气流速为2.0 m/s时，散热量为：553.95 W\n",
      "当空气流速为2.5 m/s时，散热量为：627.76 W\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEWCAYAAACEz/viAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA1r0lEQVR4nO3dd5wV1fnH8c9D772XpQiIFAFZwYLYwA62GDVWLKiJGmtiSaLGntiIGhU1Yo/YURTFAoKAgAjSe++9L9ue3x9n1lz3t7BLuXvv7n7fr9d9uXOmPXe8zDNz5sw55u6IiEjJVirRAYiISOIpGYiIiJKBiIgoGYiICEoGIiKCkoEUEjOrYmalEx3HgWZmlZP5e1mQb3xmVtPMWhRGTJKclAzkgDOzedHJ/xEzuyoqvg64Ko9lO5jZiJjpq83sljjH19rMDirAcjeZ2R35LHYu8I8C7vdkMzssV1lbMxudx7JmZg+bWZVc5RXMbNke9lHdzIbGFHUDRhUgvEbAEDPrlt+CZtbQzKoVYJtShJRJdABSvJjZh0BTYDxQF0g3swVAA2BOzHLVgeOBFKCOmZ0FfAmUBnZ7ojGzx4Dzgc25ZtUGBrj7IzHLjgOaRPPOAX4gnBzPBzLM7AMg292/jJYfBPSK2WaVUGzXx5SNdffzYqaPBc42szN3E3Kmu7eN/t4IvGhm3wA52ygL1DazRdH0Lnc/2N3dzNYAg8zs98DAaH6paPmPount7n5RzP7OBObGTJ8KjDezGtF0urvvMLPewJO5Yq0NvGtmW2PKlrv7ybmWexgYbmY/AIPcvcduvrsUIaaXzuRAMrMywCygK3AXMB94GRhLOPGnASuB3wL3ADWBo4ChwI6ovDywHZjk7n1zbf8x4GBgeq5ddwO+zCMZXAN8BjwD1ALWA3cDDwBZhGTwRLT8p8BfgAVAd3cfbmaVgaPd/UszawW85O7HRcvXBGYCNwJb8joe7j4sJp7ShKRYyd3nmlljYALQNipf7O67YpY34FPCXVVKVFwO+C8huUFINuOi5f8LnA7sjI7lMdFx3QlkADWAye5+ccw+WgAGvAWcAdxASCbzgR/cPTv2+5hZU2A04f9BE2KSgZndDsxz9w/zOhaS3HRnIAeMmfXgf1USm2Jm3Qa0BCoCw4A/uvtS4Coz6wC8Aowk3E1MI1RZfA7cn8dulhOu3pvmMW91HmV/AD4EHiecGIcD2YQTbBXg+ZhlnwWWAf+K/juccHfzN8Jdy5poWznuB94nXFG3zGPfRN8XM0sF3gP6uvvP0bx7gcfdfYuZ3QicbmYnR9OlgarARcA24FLgJkKSrAW8Ef3398C4aHsNgOPdfaKZTQYuI9xpdI9iOIVwV0Q0XQ4YDDwXE69F/70Q6Af0z/V9/gG86u5pIVf9sq07gd6EhCJFkJKBHEgTgdZ5lP8tKm8BVAdWmllD4BLClWwTwok5lXCCfR3YANSJ3YiZvUaoQlq0m/2faWbd3P26mLJNQHPCSfM1wsl/LOGK+4hc6/8cza8EXB2VbQPamtnzhKvyVGCAmR0F9AT+CgzYTTy/nHijE/QdwMdRvfx/CdVkk83sd4Qr93rAE4RnK4cAg6LYc67k7yAkgXnu3tzMntrNfnNcAdwTJZgvCFVSuwDMrCwhEWRH+8k56eec4W8GvjazK9z9P9E610Tf6caYfVSLqgbTgNPdfWc+MUmSUjKQAya6WmxOuHpeHxU3Br4HPgA6A1Xcfb2ZlSdUx9wNPODuz5nZ0YTqjXqEE3UjM3uQUDWz0N0vNbNVwKpo2/UJJ7O1hLr0Su5+Vq6wngc+JlSVOOEkdyYhSezItexdhETxZE71iLuvM7O+hLsVB+6LyseY2bHuvjHa/q+Y2RBCdVfs8flvlATTCCffHVEcu6Jt1yacwHH3aUBq9ByDaLm/AC8AHvOMYUSuXQ8xs/Qo3kOATOA7MzuckHC3RctVBQ4j3Elkm9klhGcazxLuJrLM7ObomPzHzLoQEuq/c+0vBbjX3T/IfQykaFEykHiYTbjyBrgeOAs4jnAingfg7rvMrBnwKPCpmf2WUAXxm+gquhTQnlBV8Qi/NoVw8jyUcLKbQfgt577Sh3DSGxktXxVIj5Y7iHCX8nLOgu7+BzN7FVgQWwUSY7S7vx+z/MY9HIONhJNvbuMIdxSbgFuAF4ELCHcKF7j75XltzN2fiR7sXkeokvoJeMXdh+RatG9MNVGWuy+OvtOtwFKiB+/uvsHMugLnRwk3J2YID6jfcvdrov8vuPtPUVXX3bn2N02JoHhQ01KJh22EOvdlhOoPCPX1rYHYk9cAQjXLbYQqpsOAFWZWj/B8YQKQ4u6xrVsAOhHuMhoS7jw6ExJDXvoTrmjnxZS9GLXw+Wvuhd39Mndv7u7NCQ+dB+VMxz54LYAFhOqv3E4nVH8ZoeVUQZQ1s/sID8N/Q0hoNwLXmtljtpvMFePfQBfCsfrluYq7r3X3Z4AjCc9DTiE8K1hJlIA9poWJu2eQDzO7L6pCkyJGyUDiIZXwcPRe/lfvXx+oEM3L0Rv4I+EK+gNClUNXQsLIIlwx/yln4ehK9lNC4phIeJi8JPp7PDDCzDrniuXVaJ87CQmjFbAxeth9OOEhcTyMA3qaWe+cpqlRS6vfElo3NSBUF/2KmbU0s/bR380Iia8f4S7meGAdQPQAvi+hxc8emwS6+yp3P4nQamlBHvNnEB60jyFUeZ3p7gsL8B2dmKqw6G7uLP5XFSVFiKqJJB7ecPdrAcxsDuGh73vA5cBtZtaP0Ob/BeAMd98YPVy9x93PjqomziY8YJ5oZu8DHYHHcu2nOuGE1DCmrJeZ3eDun0TTawgJpzvhLuFBwknx8mj+aWbWy91nm9nYXNuqBpQys8tjypYVsF39N4RWOusJ1WYQEtsowt3SXfy6FQ+Eu4VTgLpm1pJwZ3IX4RnKBYSkVx5obGaLCXdPFc1sobtPjLYx1MxyHkZjZg0ID47XE+4AfvkuZtaIcFd1GnAy4XlBLWCwmb1HaEE1xd1zP1vJsRwobWbrCNV1ZYFRMa2lpChxd330OWAfQrPP5wmtYJYDkwgPOW+K5jcBFhKu0jvnWrdS9N/KMWUdgHL7GMu4KI6U3W0DKB3HY3Eu4c6lHqGN/1hCIlpCeOGrFOHkP4hQhbaccOXemXAXVXcv9zcCSI3+To220SY6DmsID+pzlj0Z+BF4mtBCqELMvHqEaqh3ge+A8jHz7gWuT/TvTJ8D/9FLZ3JARVUFpdw9M6aslrtv2N10SWNmdd19baLjEImlZCAiInqALCIiSgYiIkIRbU1Up04db968eaLDEBEpUn788cd17l43r3lFMhk0b96ciRMn5r+giIj8ImqSnCdVE4mIiJKBiIgoGYiICEoGIiKCkoGIiKBkICIiKBmIiAhKBiIiRcKuzCxeG7uIYdNW5b/wPiiSL52JiJQUGVnZvP/jMp7+Zh7LN+3knMMac0qHBgd8P0oGIiJJKCvb+Xjycp76ai5LNuygc9MaPHJuR3q0qpP/yvtAyUBEJIlkZztDp67kqa/mMH/tdto3qsbLl6VyQtt65D/c9b5TMhARSQLuzpczVvPk8DnMWrWVNvWr8PzFh3FSuwaUKhW/JJBDyUBEJIHcnRFz1vLEl3OYunwzLepUZsAFnTnj0EaULoQkkEPJQEQkQcbMW8djX85m0pJNNKlZkX/+5lDO7tKYMqULv6GnkoGISCGbsGgDj385m3ELNtCgWgUePLsD53VtSrkyiWvtr2QgIlJIpizdxOPD5/DdnLXUqVKee/q048JuKVQoWzrRoSkZiIjE24wVW3hi+By+mrmampXKcuepbbn0yOZULJf4JJBDyUBEJE7mrt7KU1/NZejUlVStUIZbe7ehX48WVCmffKfe5ItIRKSIW7RuOwO+nstHk5dTqWxpbjihFVf1aEn1SmUTHdpuKRmIiBwgyzbu4Omv5/HepGWULW3079mSa3oeRK3K5RIdWr6UDERE9tOqzWk88+1c3pmwFMO49MhmXHfcQdSrWiHRoRVYXJOBmd0O9AWqAU8DXwEvR/udCfzB3bPM7DLgmmi1R9x9SDzjEhE5ENZu3cVzI+bzxg+Lyc52zj+8Kdef0IqG1SsmOrS9FrdkYGaHA8cAxwJVgD8DjxJO9sPNbCDQ18xGArcDqUA5YJyZDXP39HjFJiKyPzZuT+eF7xbw6phFpGdlc06Xxtx4Ymua1qqU6ND2WTzvDE4HZgEf8b9k8C5wQTT/M6AnsAMY6+5pQJqZzQHaAj/HMTYRkb22eWcGL49eyH9GL2R7eiZ9OzXijye2pmXdKokObb/FMxk0BNoApwDNgE+Asu7u0fxtQO3osylmvZzyXzGz/kB/gJSUlLgFLSKS2/ZdmQwas4gXRs5nS1omp3ZowM2929CmftVEh3bAxDMZpAGfufsuYI6ZbQOamJlFCaEGsB7YCMQe0ZzyX3H3gcBAgNTUVM89X0TkQNuZnsUb4xbz3Mj5bNieTq9D6nFTrzZ0aFw90aEdcPFMBmOBfmb2OFCfUFU0jFA1NJJQjTQEGA88amZlgcpAa2B2HOMSEdmjXZlZvP3DEp4dMZ+1W3dxTOs63NK7DV1SaiY6tLiJZzIYDHQFxkTT1wGLgJfMLKc10ZCoNdEAYEQUz+3R3YSISKHKyMrm3YnLeOabuazYnEa3FrV49neH0a1FrUSHFndxSwbunk1oJZTbCXks+zKhyamISKHLzMrmo8krGPD1HJZu2EmXlBr84zedOLpV7biOLpZM9NKZiJRY2dnOp1NX8tTwOSxYF4aYfOXyDhx3cN0SkwRyKBmISInj7nwxPQwxOXv1Vg6uX5XnL+7Kye3rl7gkkEPJQERKDHfn29lreGL4HKYt30LLupX514VdOKNjw0IZZziZKRmISLHn7nw/bz2PD5/NT0s20bRWRR47rxNndW6UkCEmk5GSgYgUa+MXhiEmf1i4gYbVK/DQ2R05L7UJZZUEfkXJQESKpZ+WbOSJ4XMYNXcddauW594+7bggSYaYTEZKBiJSrExbvpknh8/h61lrqFW5HHefdggXH9EsqYaYTEZKBiJSLMxZvZUnh8/h82mrqFahDLeffDCXHdU8KYeYTEY6SiJSpC1Yu40BX89lyJQVVC5XhhtPbM2VPVpQvWLyDjGZjJQMRKRIWrphB//6ei4f/LSccqVLcU3Pg7imZ0tqFoEhJpORkoGIFCkrN+/kmW/mMXjiUsyMy45sznXHHUTdquUTHVqRpmQgIkXCmq1pPDdiPm/+sAT3aIjJ41vToHrRGWc4mSkZiEhS27A9nRdGzufVsYvIyHJ+c1gTrj+hVZEeYjIZKRmISFLavDODl0Yt4D+jF7IjI4uzOodxhlvUqZzo0IolJQMRSSrbdmXyyuiFvDhqAVvSMjm9Y0Nu6tWa1sVoiMlkpGQgIklhZ3oWr41dxPMj57NxRwa9DqnPzb1b075R8RtiMhkpGYhIQqVlZPH2+CU8++181m3bRc82dbmldxs6N62R6NBKFCUDEUmI9Mxs3v1xKc98M4+Vm9M4omUtnrv4MA5vXvyHmExGSgYiUqgys7L54Kfl/OvruSzbuJPDUmrw+HmdOKpVnUSHVqIpGYhIoXB3Pvl5JU8On8PCddvp2Lg695/VgePalLwhJpNRXJOBmU0AtkeTC4EPgAHAkqjsHncfaWZ3AH2islvdfVw84xKRwjVt+WbuHTKdiYs30rZBVV64pCsntSu5Q0wmo7glAzMrB5Rx9+Niyu4F7nD3wTFl7YDTgB5ACvA+kBqvuESk8KzftovHvpzNfycspValcjx6bkfO69q0xA8xmYzieWfQEahsZsOBUsDdQFegu5ndAIwH/gwcA3zh7g4sNrMyZlbN3bfEMTYRiaOMrGxeH7uYJ7+aw870LK44ugU3nthaPYkmsXgmgzTgcWAg0BoYBvwb+BBYADwP/B6oAmyKWW8bUBv4VTIws/5Af4CUlJQ4hi0i+2P03HXc98l05q7ZxjGt63BPn3a0qqcXxpJdPJPBXGBedMU/x8zWAYPdfQmAmX0MnAP8CMS+VVIDWJ97Y+4+kJBYSE1N9TjGLSL7YOmGHTwwdAZfTF9NSq1KDLykK731XKDIiGcyuJxQLXSNmTUGqgHfm9nRUULoTagqGgP8C3jEzFoAGaoiEik6dqRn8tyI+bzw3QJKm3H7yQdzZY8WGmu4iIlnMhgEHGNmowEH+hGqhN41s13AdOBld88ys2/NbBRQGrgujjGJyAGS01T04c9msnJzGmd2bsQdp7alYfWKiQ5N9kHckoG7pwOX5DFreB7L3g/cH69YROTAmr5iM/cNmcH4RRto36ga/7qwi94cLuL00pmIFNiG7ek8/uVs3h6/hBqVyvHwOR35bWpTSqupaJGnZCAi+crMyubNH5bw+Jez2Z6exWVHNeemE9tQvZKaihYXSgYiskdj5q3jvk9mMHv1Vo5uVZt7+rSnjcYWKHaUDEQkT0s37OChz2by+bRVNKlZkecv7srJ7dVUtLhSMhCRX9mZnsVzI+fzwsj5mMGtvdtwdc+WaipazCkZiAgQmop+NnUVDw6dwYrNafTp1Ig7T21LoxpqKloSKBmICDNXbuHeIdP5YeEGDmlYjSfP70z3lrUTHZYUIiUDkRJs4/Z0nhg+hzd/WEy1imV54KwOXNgtRU1FSyAlA5ESKDMrm7fHL+Hx4XPYsjODS45oxs2921CjUrlEhyYJomQgUsKMW7Cee4dMZ9aqrRzZsjb39G1H2wbVEh2WJJiSgUgJsXzTTh76bCZDf15J4xoV+fdFh3FqhwZqKiqAkoFIsZeWkcULIxfw3Mh5uMNNvVpzTc+DqFhOTUXlf5QMRIopd2fYtFU8MHQmyzft5PSODbnztLY0qVkp0aFJElIyECmGZq/ayn2fTGfM/PW0bVCVt68+giMPUlNR2T0lA5FiZPOODJ78ag6vj1tMlfJl+PuZ7fldtxTKlC6V6NAkySkZiBQDWdnOfycs4bEvZrN5Zwa/657Crb0PpmZlNRWVglEyECniJizawD0fT2fGyi10a1GLe/u0p10jNRWVvaNkIFJErdy8k4c/m8WQKStoVL0CT1/YhTMObaimorJPlAxEipi0jCxeGrWAZ7+dT5Y7N57YmuuOVVNR2T9xTQZmNgHYHk0uBO4DXo72OxP4g7tnmdllwDXRco+4+5B4xiVSFLk7X85YzQNDZ7B0w05Oad+Au08/hKa11FRU9l/ckoGZlQPKuPtxMWXvEE72w81sINDXzEYCtwOpQDlgnJkNc/f0eMUmUtTMXb2V+z6Zweh562hTvwpvXtWdo1vVSXRYUozE886gI1DZzIYDpYC7ge7ABdH8z4CewA5grLunAWlmNgdoC/wcx9hEioTNOzMY8NVcXh27iMrlSnNPn3ZcfEQzyqqpqBxg8UwGacDjwECgNTAMKOvuHs3fBtSOPpti1ssp/xUz6w/0B0hJSYlb0CLJICvbGTxxKY99MZsNO9K54PAUbjupDbWrlE90aFJMxTMZzAXmRSf/OWa2DjjMzCwqqwGsBzYCsaNr55T/irsPJCQWUlNTPfd8keJi4qIN3PvJdKYt38LhzWvyap9udGhcPdFhSTEXz2RwOdAVuMbMGgPVgCGEqqGRwOnR9HjgUTMrC1Qm3EXMjmNcIklp1eY0Hvl8Jh9NXkGDahUYcEFn+nZqpKaiUijimQwGAceY2WjAgX7AauAlM8tpTTQkak00ABgRxXO7u++KY1wiSWVXZhYvjVrIs9/OIzPLuf74Vlx33EFULq+W31J4CvRrM7MKwBnAUUBNwkn9W2C4u2fntU7UGuiSPGadkMeyLxOanIqUGO7O1zPXcP/QGSxev4Pe7erzl9MPoVntyokOTUqgPSYDC/enNwLnEKp2vic87K0FHA/cZmZP670Akb0zb802/v7pDL6bs5aD6lbmtSu60bNN3USHJSVYfncG7YBZ7n5sHvPeNbNSwGVmVtXdtx748ESKly1pGfzrq7kMGrOIimVL89cz2nHpkWoqKom3x2Tg7tOB6XuYnw28cqCDEilusrOd935cxj++mMX67emcn9qU204+mDpqKipJQk+oROJs0pKN3DdkOlOWbeawlBq8cnk3OjZRU1FJLkoGInGyZksajwybxQeTllOvanmePL8TZ3VurKaikpTyTQZmdqi7q2sIkQLalZnFK98v4umv55KR5Vx33EH84fhWVFFTUUliBfl1DjCzFEJromHAl+6+Ka5RiRRR38xazd8/mcGi9TvodUg9/nJ6O5rXUVNRSX75JgN3P97MygNHAscBV0etiEa6+9/jHJ9IkbBg7Tbu/3QG385eS8u6lRnU73COO7heosMSKbAC3be6+y4z+5HwfkFVoAvQOY5xiRQJW9MyeOabefzn+4WUL1Oau087hMuOak65MmoqKkVLQZ4Z/Ak4lZAEvgI+Be5w94w4xyaStLKznQ9+Ws6jw2axdusuzuvahNtPOZh6VSskOjSRfVKQO4O7gM+BhwhVQxp0Rkq0yUs3ce+Q6UxeuonOTWvw4qWpdG5aI9FhieyXgiSDeoRnBX2BJ81sAWFgms/cfUkcYxNJKpt3ZvDg0BkMnriMulXL89h5nTinS2NKlVJTUSn6CvIAOR34EvgyenB8CXAn8CygEbilRBg9dx23vTuFtdt20b9nS244oRVVK5RNdFgiB0xBnhkcDhwbfQ4GJgCPEBKESLG2Mz2LR4fNYtCYRRxUtzIDLz2KQ5vUSHRYIgdcQaqJHiac+P/m7j/FOR6RpDFl6SZuHjyZBWu3c/lRzbnj1LZUKKubYSmeClJN1KswAhFJFhlZ2TzzzTye+XYe9aqW540ru9OjdZ1EhyUSV3o/XiTGvDXbuGXwZH5etpmzuzTm3r7tqV5Rzwak+FMyECG8N/Dq2EU88vksKpYrzb8vOozTOjZMdFgihSa/kc7K5zcesZmVBnx3w1+KJLsVm3Zy+3tT+H7eeo4/uC6Pnnso9arp5TEpWfK7M2hqZg8AbwND3T0zZ4aZVQPOJwx/eQWQFrcoReLA3flo8nL+9vF0srKdh87uyIXdmqqLaSmR9tiBirvPA64kDH85ycx+MrPvzGw68AWQDVzi7rtNBGZWycwWmFlbM+toZivMbET0OT9a5jIzGxN9+h64ryeStw3b0/nDW5O4+Z0ptKlflc//eAy/656iRCAlVkFaE20nNC992MxqAjWANVF5QdwfrQPQFXjK3f+RM9PMagG3A6lAOWCcmQ1TtxcSL9/MWs2f35/Kph3p/OmUg7mm50GU1lvEUsLt1QNkd98IbCzo8tELa7WAnMFxugIdzOwMYC5wE3A4MDa6u0gzszlA25h1RA6I7bsyeWDoTN4ev4SD61dlUL/Dad9Iw0+KQBxbE5lZGcKbyucD70XFE4HX3H2Cmd0N3BeVbYpZdRtQO4/t9Qf6A6SkpMQrbCmmJi7awC2Dp7B04w6u6dmSW05qQ/kyeoFMJEc8O13/E/C6u6+LKfvQ3Sfk/A0cSrjTqBqzTA1gfe6NuftAd09199S6devGKWQpbnZlhu4kfvvCWLLdeaf/kdx52iFKBCK5FPjOwMyqu/vmmOkr3P0/e1jlFCDbzC4nDITzGlA1Wm8s0BsYH30eNbOyQGWgNTB7b7+ISG6zVm3hpv9OZtaqrZyf2pS/9mmncYhFdmNv/mW8Z2YPR383Bq4DdpsM3L1nzt9mNgK4FqgIPG1mmcAq4Gp332pmA4ARUTy35/dug8ieZGU7L41awONfzqFaxTK8dGkqvdrVT3RYIkktv5fOlgEO/BQt+y+gEqFeP3MPq/6Kux8XM9kjj/kvAy8XdHsiu7N0ww5uHTyF8Ys2cHL7+jx0dkdqVymf6LBEkl5+dwYLCSf99YR6/ZlAnegjkjTcncETl/L3T2ZQyozHz+vEOYc11nsDIgVUkGoijz4ALYAq0ScrXkGJ7I21W3dx5wc/89XMNRzZsjb/PO9QmtSslOiwRIqUgj4zyLm8yiK8dZzzEUmoYdNWcdeHU9m2K5O/ntGOfkc11zCUIvugoMkg585gMaGKqC6QEZeIRApgS1oG9w2ZwfuTltGhcTWe/G1nWtevmv+KIpKn/JJBM0Ii2BpNH0p4gLx5t2uIxNmY+eu4/d2fWbUljRtPaMX1J7SmXJl4vjIjUvztMRm4+y+v+prZcOCiaLIpcHcc4xL5f9IysvjnF7N5efRCWtSpzHvXHkmXlJqJDkukWNib9wzOdPcd0d8/mlm5eAQkkpdpyzdz8zuTmbtmG5ce2Yw7Tm1LpXJ6gUzkQCnwvyZ332FmxwML3X0RMMbMurn7+LhFJyVeZlY2z42Yz4Cv51K7Sjleu6IbPduoOxKRAy3fZGBmjQhvDtcALgdam9ldwF+BW+IZnJRsC9Zu45bBU5i8dBN9OzXi/jM7UL2SxiMWiYeC3BmMAb4j9Cw6hdDtdC/gdNSiSOLA3Xlj3GIe/Gwm5cuU5ukLu9CnU6NEhyVSrOXXHcVRwA7gJaAacL27P2FmPwCfAQ8A38Q9SikxVm1O4/b3pjBq7jp6tqnLP39zKPU1HrFI3OV3Z3ASoSnp8cANwNyo2igVuNLdF8Y5PilBhkxZwV8/mkZ6Zjb3n9WBizUMpUihyS8ZPEjoVnoooV+i1cAHwLHA/dHwlLozkP2yaUc6f/loGp/+vJIuKTV44redaVGncqLDEilR8ksG7QjvFHQADgJGA1cBZwDjlAhkf42YvYY/vfczG7anc9tJbbj22IMoU1ovkIkUtvySwb2Et44dOIGQEKoCLxIGpRHZJzvSM3nos5m8MW4JretV4T+XH06HxhqPWCRR8rsEuw74zt0vBN4Cnid0a70K+MjMeu5pZZG8TFqykdP/NZo3f1jCVT1a8MkNPZQIRBIsvzuDjcBzZlYbeA7Y6O43ApjZqcD1Zjba3dWDqeQrPTObp7+Zy7PfzqNh9Yq8ddURHHlQ7USHJSLk3zfRLmC8mS0F/gnMjIawLAc8RujSuhTqzlryMXf1Vm4ePJlpy7fwm65N+FufdlSroBfIRJJFfu8ZlHb3LEIvpfOAw4CBhGTwJ3d/M/4hSlGWne385/uF/OOL2VQpX4bnL+7KKR0aJDosEcllt8nAzCoC75hZdaA+oZfSbYTmpuuB28xsvruPK5RIpchZtnEHt707hXELNtDrkHo8fM6h1K2q8YhFktFuk4G77wT6mlkq4V2DR4AlwOBokZnAv83s9GjZPJlZJWAacBqwBhhEaJG0Fujn7tvN7CTgPkJ102vu/sL+fjFJHHfn/UnLuW/IdLLd+ce5h3JeahO9QCaSxPKrJrob6A98DnwK3EpoYlqD0GdR9T0lgsj90fIAfwY+cfcXo87u+pvZ08AAoAewhfCM4kN3X7NP30gSav22Xdz14VS+mL6abi1q8fh5nWhaS+MRiyS7/JqWjiNUDc2IpmcT+iSa7e59CB3X7ZaZHQ7UAn6OinpG6xP99xjCuwvL3H29u2cAI4Ej9vJ7SBIYPmM1Jz/1Hd/OWstdp7Xl7auPUCIQKSLya1rai9BBXSohEVQEagIVzawdoUVRnsysDKFq6Xzgvai4NqH3UwhJpnaustjy3NvrT7hLISUlJfdsSaCtaRnc/+kMBk9cxiENq/HmVZ05uIHGIxYpSvJLBh8Rup4YBdQDhkflK4F+gJvZHbt5z+BPwOvuvi6mrngj4XnBdkLV0fqYshw55b/i7gMJLZlITU31fOKWQvLDgvXc+u4UVmzaye+PO4iberXReMQiRVB+7xn8YGadCFU99d19+l5s+xQg28wuBzoDrwFTgZOBVwnjIYwC5gLNzKwGobvsnoTnDJLE0jKyeGL4HF4ctYCUWpV499oj6dqsVqLDEpF9lO/gNtFV/7roU2Du/ktXFdGLatdG23jVzPoRWhb1c/dMM7uV8AyhFPCsu6/em31J4Zq+YjO3vDOF2au38rvuKdx92iFULq/xiEWKskL5F+zux8VMnp7H/M/434NlSVJZ2c7zI+fz1FdzqFGpHK/0O5zjD66X6LBE5ADQ5ZwUyKJ127n13Sn8uHgjp3dsyANndaBm5XKJDktEDhAlA9kjd+et8Ut4cOhMypQyBlzQmb6dGukFMpFiRslAdmvNljT+/P7PfDt7LT1a1eGf5x1Kw+oVEx2WiMSBkoHk6bOpK7n7w6nszMjivr7tueSIZpQqpbsBkeJKyUB+ZfPODO75eBofTV5BpybVeeL8zhxUt0qiwxKROFMykF+MnruO29+bwpqtu7i5Vxv+cLzGIxYpKZQMhJ3pWTw6bBaDxizioLqV+fD3R3FokxqJDktECpGSQQk3Zekmbh48mQVrt9Pv6Ob8+ZS2VChbOtFhiUghUzIooTKysnnmm3k88+086lUtz5tXdefoVnUSHZaIJIiSQQk0b802bhk8mZ+XbeacLo25p297qlfUeMQiJZmSQQni7rw6ZhEPfz6LSuVK89xFh3Fqx4aJDktEkoCSQQmxKzOLO96fyoc/LeeEtvV45NyO1KtaIdFhiUiSUDIoATZsT+ea1ycyYdFGbu3dhutPaKXuJETkV5QMirl5a7Zx5asTWLk5jacv7EKfTo0SHZKIJCElg2JszLx1XPvGj5QtXYq3rz6Crs1qJjokEUlSSgbF1DsTlnD3h9NoUacy/7n8cA1MLyJ7pGRQzGRnO49+MYsXRi7gmNZ1ePaiw6hWQc1GRWTPlAyKkZ3pWdz8zmSGTV/FxUekcG+f9upbSEQKRMmgmFizJY2rXpvI1OWb+dsZ7eh3dHO1GBKRAlMyKAZmrNjCla9OYPPODF68JJVe7eonOiQRKWLiVodgZqXM7EUzG21mo8yso5n1MbMFZjYi+hwbLXuHmX0ffY6IV0zF0TezVnPe82MAePfaI5UIRGSfxPPOoA9Q2t17mNnxwEPAj8Ad7j44ZyEzawecBvQAUoD3gdQ4xlUsuDuDxizi/k9n0L5RdV66LJX61fRGsYjsm7glA3f/2MyGRpPNgfFAN6C7md0QTf8ZOAb4wt0dWGxmZcysmrtviVdsRV1mVjb3fTKD18ct5qR29Xnqgs5UKqcaPxHZd3FtauLumWb2CvAUsB74Grge6AlUAX4P1AY2xay2LSr7FTPrb2YTzWzi2rVr4xl2UtualsGVr07k9XGLuebYljx/cVclAhHZb3E/i7h7PzO7A5gEtHP3zQBm9jFwDqHqqHrMKjUIiSP3dgYCAwFSU1M9zmEnpWUbd3DloInMX7uNR87pyAXdUhIdkogUE/F8gHxplAQAdgBbgRlmlnMG602oKhoFnBSt0wLIUBXR//fTko2c9ez3rNi8k1ev6KZEICIHVDzvDN4HXjOz7wAD/ghkA++a2S5gOvCyu2eZ2bdmNgooDVwXx5iKpE9/XsGtg6dQv1oF/tv/cFrVq5LokESkmInnA+TtwLl5zBqex7L3A/fHK5aiyt159tt5PPblHFKb1WTgpanUqlwu0WGJSDGkJ49JaldmFnd+MJUPJi3nrM6NePQ3h1K+jAaqF5H4UDJIQhu3p3PN6z8yftEGbu7VhhtP1GA0IhJfSgZJZv7abVw5aAIrNqcx4ILOnNm5caJDEpESQMkgiYydv55r3/iRMqWMt6/uTtdmtRIdkoiUEEoGSWLwxKXc9cFUmtepzCsajEZECpmSQYJlZzv//HI2z42YzzGt6/DM7w6jekUNRiMihUvJIIF2pmdxy+DJfD5tFb/rnsJ9fdtTVoPRiEgCKBkkyJotaVz92kR+Xr6Zv5x+CFf2aKEWQyKSMEoGCTBz5RauHDSBTTszGHhJKr01BoGIJJiSQSH7dtYarn9rElUrlGXwNUfSoXH1/FcSEYkzJYNCNOj7hfz90xkc0rAaL192OA2qazAaEUkOSgaFIDMrm/s/ncGrYxfTu119BmgwGhFJMjojxdnWtAxuePsnRsxeS/+eLfnzKW0pXUoPikUkuSgZxNHyTTu5ctAE5q7ZxkNnd+R33TUGgYgkJyWDOJm8dBNXvTqRXZlZvNqvGz1a10l0SCIiu6VkEAefTV3Jze9Mpl618vy3f3da1aua6JBERPZIyeAAcnf+PWI+//xiNl2b1WTgJV2pXaV8osMSEcmXksEBkp6ZzZ0fTOX9Scs4s3MjHj33UCqU1WA0IlI0KBkcABu3p3PtGz/yw8IN3NSrNX88sbW6lhCRIkXJYD8tXLedKwZNYPnGnTx1fmfO6qLBaESk6IlbMjCzUsALwCGAA78HtgIvR/udCfzB3bPM7DLgmmjVR9x9SLziOpDGLQiD0ZQy462ru5PaXIPRiEjRFM87gz5AaXfvYWbHAw8BOwgn++FmNhDoa2YjgduBVKAcMM7Mhrl7ehxj22/vTlzKXR9OJaVWJV65vBsptTUYjYgUXXFLBu7+sZkNjSabA+OBK4ELorLPgJ6EBDHW3dOANDObA7QFfo5XbPsjO9t5fPhsnv12Pke3qs2/L+qqwWhEpMiL60gq7p5pZq8ATwHrgbLu7tHsbUDt6LMpZrWc8l8xs/5mNtHMJq5duzaeYe9WWkYW1789iWe/nc+F3VIY1K+bEoGIFAtxf4Ds7v3M7A5gElDJzCxKCDUICWIjEPtWVk557u0MBAYCpKameu758bZmaxpXv/YjPy/bpMFoRKTYidudgZldGiUBCFVBW4FvCFVDAKcDowjVR0eZWVkzqwG0BmbHK659MWvVFs5+dgxzVm3lhYu7ctUxLZUIRKRYieedwfvAa2b2HWDAH4G5wEtmltOaaEjUmmgAMCKK53Z33xXHuPbKiNlruP6tn6hcvjTvXqvBaESkeIrnA+TtwLl5zDohj2VfJjQ5TSqvjV3EvUOmazAaESn29NJZHrKynfs/ncGgMYvodUgYjKZyeR0qESm+dIbLZduuTG58+ye+mbWGq3q04M7TDtFgNCJS7CkZxIgdjObBsztwUfdmiQ5JRKRQKBlEpizdxFWvTSQtPYtB/Q7nmNZ1Ex2SiEihUTIAhk1byU3vTKZOlfK8dVV3WtfXYDQiUrKU6GTg7jw3cj7/GDabw1JqMPDSVOpoMBoRKYFKbDJIz8zm7g+n8u6Py+jTqRH//I0GoxGRkqtEJoNNO8JgNOMWbODGE1tzcy8NRiMiJVuJSwYL123nykETWLZxJ0+e34mzuzRJdEgiIglXopLBlKWbuOyV8Rjw5tXdOVyD0YiIACUsGTSpWZHOTWtwX9/2NKtdOdHhiIgkjRKVDGpXKc+gft0SHYaISNKJ6+A2IiJSNCgZiIiIkoGIiCgZiIgISgYiIoKSgYiIoGQgIiIoGYiICGDunugY9pqZrQUWJzqOPagDrEt0EHug+PaP4ts/im//7E98zdw9z5G7imQySHZmNtHdUxMdx+4ovv2j+PaP4ts/8YpP1UQiIqJkICIiSgbxMjDRAeRD8e0fxbd/FN/+iUt8emYgIiK6MxARESWDvWZmd5jZ99HniJjyxmY2Iuazwcz+GM2bEFP+SiHE2NPMvsujvIuZjYo+9+X3nRIQ3/lm9oOZjTGz582sVFS+NOb4PZzA+G4wsxkxsRxswQAz+87MvjWz1omILxl+f2ZWxswGRb+v8WbWJ9f8k8xsbPQ7uyYqK9TjV4AY/xj9Bseb2d+ismpmti7mGP4xgfE9bmaTYmKpbmYVzOz1aPpzM6u3Tzt3d30K+AHaAd8BBjQDJu5mua7RcmWBcsBPhRjjn4GfgXF5zBsHtIn+/gLoUtDvFO/4gArAQqByNP0O0AdoDnyWJMdvENAtV9mpwJvR30cBHycqvkT//oDLgGeiv+sAi2PmlQFmArWjuH4C6iXg+O0pxuZRXKUJF8pjgY7AccC/E30Mo7IRQL1cZdcBD0Z//w4YsC/71p3B3jkG+MKDxUAZM6uWx3IvADe4ewbhx1TZzIab2deFcOU9Dzg3d6GZlQdqufucqOhzwvcp6HeKa3xAOnCUu2/PCRnIJJzYGpnZN2Y21MwOjmNse4qPKJY7zGy0md0ZlfUEPgNw9zHAoQmML0eifn/vA3dHf2fnmncQsMzd10dxjQSOoPCP355iXA6c7O5Z7p5NSAg5v8EuZjbSzAabWcNExGdmBhwMvBD9BvtFs345htF/j9mXHSsZ7J3awKaY6W1R2S/M7AxggbtPiYrSgMeBkwgZ/C0zi9two+7+PpCRx6zawOaY6ZzY8/1OhRGfu2e7+0oAM7sRqAEMA1YDD7v7CcDDwJvxim1P8UXeBa4BTgB6mFlf/v/xK5NTvZWA+BL6+3P3be6+2cyqAu8Bf42ZvbvfWWEfv93G6O4Z7r4mqrp6DJjm7jOBOcA97n4s8DHwbCLiAyoB/wYuAk4B/mBmnfj1Mdznf78lagzkA2AjUD1mugawPtcylwCvxkzPBeZ5uIebY2brgAbAsjjGmZeNQNWY6RqE2DPI/zsViujK52GgPXC2u7uZ/Qj8AODuo82soZlZdDwLO7Yn3X1rND2UcBWb+7h6dFWZKAn9/ZlZY+AD4EV3fy1m1u5+f4V+/PYQI2ZWDngZ2ElI/ADfEJIqwIfAfcTRHuJLA55w9x3Rct8Q7vxij2EN9vHfr+4M9s4owhUWZtYCyHD3LbmWOR4YHjN9OfCvaJ3GQDVgRdwjzcXddwKbzaxldGI7FRhNwb5TYXmBcHzOjKku+itwWxRfF2BRYSeCSBVgVvQw0YBewHjC8Ts5iu9YQp1zIiXs9xdVn3wJ3O3uL+WaPRdoZmY1ohNuT0KSL9Tjt6cYzaw04cp/urv3d/fMaNaLwG+jv3sT/r8XenxAK+CH6CFzWUJ10I/EHEPg9Gh6r+nOYC+4+7SoxcMowkOm68zsEqCcu79sZrWADe6eHrPaIOAYMxsNONCvMK8cY+MDrgdeI9THf+Xuk6JlfvWdCiu22PiAScBVhB/yN+F8ywDgH8AbZjaSUH97RSLii/7//gn4mvB8Y7i7fxlVaZwaxQdwdQLjS/Tv7w5CFcVfzOwvUdnXwCR3H2pmtxLqtEsBz7r7ajP7jMI9fruNkfA7PB6oaGanRPPuBO4C/mOhBdS2OMeY3zF8hfBgOxN4xd1nmtki4BUzG0G4e7hsX3asl85ERETVRCIiomQgIiIoGYiICEoGIiKCkoGUEGYWtxeFRIoDJQMpNszsWjO7II/yqsAlZjZuN5/q0Vun15tZ/Zj1zjKz23Jtq6KZPW1mlcysgZk1t9AB25vR381394avmbU2s4oH/Iv/eh/x7s5Biik1LZUizczWAvMJbcTTgcqE92c2Ezr0amlmFxP6wflvHpuY6u6bo22lAvcS2uPXBToDk2OWPdvdV5rZSYT+apYROverF21/SLTc0znbjImzCvCyu5+/j9/zHcI7AjvyWe5B4EN3n7gv+5GSSy+dSVE3kfAi2iB3P9nMfgM0cPdnzOyT6Cr9ZsLLa2fksf4KwpvZRxJe9jnT3bPM7EJCdxMPAX8DHnH39dFdxhJgFtCI0H3GekIXD+cAc3IngsiN5J2M8mVmFYBS+SWCyLOEN7n75LegSCxVE0lRdz2hQ6+ctzV3EXqYPINw5f4HQnfYpxO6Io79bHf3BdF6kwn9+pxsZo8SuiDoQegj5nDg3agbiu6Et3pz+owZROg47EOg7x7ivBAYmrvQzC43s/cs9EM/w8yuMLOPzWxelNggdIz3rZnVNbOvLIz3MMHMOufenruvAKrHVneJFITuDKTIivpn+StwpbvPioq/IlzhNwE+InSVnJ67CwYLXWHflTPt7jvN7DqgE6EjsucIXQ9kEKqMOkV9In0FfBW9+r+V0EfRncA6QrVSbL9AOfuqBWTn6iYiVjV3Pyk6+f+JkHCOif5+j5DIHiVURa0nXPW3JfQzlJefgSOj7y9SIEoGUmS5e0bUf8t7UV9G9Ql9wK+NFunl7mm7WX0doYfH2O1tMLOLgAnA88CnwCrgenc/Lo/9TzezT4A1hOcGq4DX89hXDULi2J2c7qY3ATOi3lo3EAb8AUhx9yVmlvOM4gPCc40Hd7O9TUDNPexP5P9RMpAizd2XEa6YMbPrgW3uPqgA6603s1oWBvJJiTohbEroCrgc4WT7i6iKiJgeU83MBhDuHP4C/B54AHjNzM7O9dxgDVBrD+Fk7W5G1DpoajR5HLDe3U81s6OA+4ET81itFjBmD/sT+X+UDKTIirr7/TCmqD6QbWbXxpQd7e67O9nOAa4kJIBpwD2EZwC/jebFOgTob2afE07+64BjCdU25xKGTDyM0NV1T+CTnBXdfZuZbTWzGu6+aS+/5hmEOxQI1T9vRr1nVgT+vpt1DiVUXYkUmJqWSollZi0Jdfx9CNVLfyOME30B4WHwBYSqn56EK/C/ELoTzooZb6Gg+7qa8O9t4AH7AnnvpxWhL/x++S4sEkPJQCRG9D7A9gM9gE7UxPUN4JJoDOC4MLNngH+4+5J47UOKJyUDkUJiZpUILZsy81143/dRfTfvOYjskZKBiIjopTMREVEyEBERlAxERAQlAxERQclARERQMhAREeD/AMeGq/YlGmaVAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "d = 7.5e-3\n",
    "distance = 10e-3\n",
    "t_f = 25\n",
    "t_w = 120\n",
    "W = 10e-2\n",
    "\n",
    "L = 6e-2  # 圆柱的高度。缺少条件，这里自己取的6 cm。\n",
    "lambda_ = 180  # 导热系数。缺少条件，这里自己取的180 W/(m-K)。\n",
    "\n",
    "\n",
    "def h(u):\n",
    "    return 5.12 * u**0.65 / d**0.35\n",
    "\n",
    "\n",
    "u = np.array([0.6, 1.0, 1.5, 2.0, 2.5])\n",
    "h = h(u)\n",
    "\n",
    "perimeter = np.pi * d\n",
    "A_c = np.pi * (d/2)**2\n",
    "H_p = L + A_c / perimeter\n",
    "eta_fin = fin_tip_efficiency(H_p, perimeter, A_c, lambda_, h)\n",
    "# print(f'肋片的效率为: {eta_fin:.2%}')\n",
    "number = 10*10\n",
    "A_f = number * H_p * perimeter\n",
    "A_r = W*W - number*A_c\n",
    "eta_overall = overall_fin_surface_efficiency(A_f, A_r, eta_fin)\n",
    "# print(f'肋面的总效率为: {eta_overall:.2%}')\n",
    "R_overall = 1 / (h * (A_f + A_r) * eta_overall)\n",
    "Phi = (t_w - t_f) / R_overall\n",
    "[print(f'当空气流速为{speed} m/s时，散热量为：{phi:.2f} W') for (speed, phi) in zip(u, Phi)]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(u, Phi)\n",
    "ax.set_xlabel('空气流速（m/s）')\n",
    "ax.set_ylabel('散热量（W）')\n",
    "ax.set_title('散热量随空气流速的变化')\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 习题02-93"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEFCAYAAADjUZCuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAkLUlEQVR4nO3dd3xUdb7/8dcnkwZJaEnohICg0gWCQIBYFxHbtaACUhRpYlt3r8vede91m1stizRRFEWxAWtDRVR6ERJAQEB6b4HQWwh8f39kvL9cFhAkM2fK+/l4zOORmTPMvI8DvnPOnO/3a845REREYrwOICIioUGFICIigApBRET8VAgiIgKoEERExC/W6wA/VVpamsvMzPQ6hohIWMnLy9vtnEs/07awLYTMzExyc3O9jiEiElbMbOPZtumUkYiIACoEERHxUyGIiAigQhARET8VgoiIACoEERHxC2ohmFmOmc04w+PNzWym//a7YGYSEZFiQRuHYGa/AroDR86weQTQ0zm3yswmm1lz59yiQORYtfMgn3y7jYQ4H/G+GBLiYkiKjyUtJYG05HjSkxNIS04gJsYC8fYiIiErmAPT1gB3AmNLPmhmCUAl59wq/0OfAR2AfysEM+sH9APIyMj4SSFW7TzIkK/XnPM5CbEx1ElLIjM1ibrpSTSpUZ4mNctTo0IZzFQUIhKZglYIzrkJZpZ5hk2pwP4S9w8Btc7yGqOAUQBZWVk/aWWfm5tW56Ym1Sg8eYrColMcLzrFwWNF7Dl0nN2HjpN/8DibCo6wfvdhVu86yJcrdlJ0qvitUpPiaZ5RkXb1UmlfL416lZNVECISMUJh6oq9QEqJ+xWAPYF8QzMjIdZHQqyPFCAtOYE6aUlnfO7xopOs3H6QJVv28e2W/SzYUMCXK3YCUDklgZxL07mxcVXa108jIdYXyNgiIgHleSE4546a2X4zqwusB24E/tPjWP8rIdZHs1oVaFarAj38j20uOMKctbuZuXo3k7/bwfi8LSQnxHLt5ZW5uWk1rrm8MnE+XcAlIuHFs0Iwsx5AvHNuNPAw8AZgwJfOuYVe5ToftSqV5Z5KGdzTKoPColPMWbubz5ft4IvlO/no222kJcdzR4uadGlZk/pVUn78BUVEQoA595NOxXsuKyvLhdpsp0UnTzF9VT7v5W7mqxW7KDrlaJVZkQfa1eFnDasQq6MGEfGYmeU557LOuE2FEBi7Dx3nXwu38sa8DWwuOEqNCmXonZ3JPVfWolxinNfxRCRKqRA8dPKU48sVOxk9az3z1xeQnBBL7+xM+rSvQ8WkeK/jiUiUUSGEiGVb9zN82ho+XbqDpHgfvbIzebBDXSqpGEQkSFQIIeb7HQcZ8vVqPl26nTJxPvrl1KVfTl3Kxnt+0ZeIRDgVQohavfMgz01ZxWfLdlA5JYEnfnYpXbJq4dO0GSISIOcqBF324qH6VVIYcV9LJgxsS82KZRg8cSk3/nMG01flex1NRKKQCiEEtKxdiQkDsxnRvQWFRafo9ep8BozNY9u+o15HE5EookIIEWbGjU2qMfnnOfznDZcxbdUurnt2OiOmraWw6JTX8UQkCqgQQkxCrI9B19Rjys+von39NP76+Uo6D5nJN+sCOr2TiIgKIVTVqlSWl3tm8WrvLI4XneSeUfP47w+Xcfh4kdfRRCRCqRBC3LWXV2Hy4zk80K4OY+dt5IYXZjB7zW6vY4lIBFIhhIGy8bH89y0Neb9/W+J9MXR/5Rt+PXEpB4+d8DqaiEQQFUIYycqsxKePdaB/Tl3eXbCJTi/MZP76Aq9jiUiEUCGEmcQ4H7/u3IAJA7OJ8xn3jprLPyZ/z4mTuhJJRC6OCiFMNc+oyKRHO3BXy5oMnbqGu0bMYf3uw17HEpEwpkIIY0kJsfztrmYM796CDXuOcNOQmby3YDPhOh2JiHhLhRABOjepxuePd6BZzQo8OWEJj7+7WJenisgFUyFEiGrly/DWg635ZcdL+fjbbdw6dBbf7zjodSwRCSMqhAgSE2M8fG193nywNfuPFnHbsFmMz9vidSwRCRMqhAiUfUkanz7WnitqVeCX73/Lr8Yv4diJk17HEpEQp0KIUJVTEnmzT2sevqYe7+Zu5j+GzWbjHl2FJCJnp0KIYLG+GH55w2W8dn8rtu8/xq1DZzNDay2IyFmoEKLANZdV5qOH21GtfCK9X5vPS9PX6tJUEfk3KoQoUTs1iQkDs+nUuCp//mwlj76zmCOFujRVRP4/FUIUSUqIZVi3FjzZ6TI+WbKNO0fMZXPBEa9jiUiIUCFEGTPjoavr8WrvVmzde4Rbhs5ijqbTFhFUCFGr+HuF9qQnJ9Dz1fm8M3+T15FExGMqhCiWmZbEhIeyya6XxuCJS/nTpOWcPKUvm0WiVVALwcwGm9ls/63Nadu6mtkiM5tlZoOCmSualUuM49VeWfTOzuTlmevpPzZX8yCJRKmgFYKZNQQ6A+2BbsDQEttSgb8A1wFXAV3NrGWwskW7WF8MT9/aiN/f1oip3+dz18i5bN131OtYIhJkwTxC6ABMdsU2ArFmVs6/rS6w2DlX4Jw7Ccz2P///MLN+ZpZrZrn5+RpgVdp6ts3k1d6t2FJwhNuGzmbRpr1eRxKRIApmIaQC+0rcP+R/DGAN0NDMqphZGYqPFMqc/gLOuVHOuSznXFZ6enqg80alqy5NZ+JD2ZSJj+HeUfP4ZMk2ryOJSJAEsxD2Aikl7lcA9gA45/YCjwHjgXeA+cD2IGaTEupXSeGDh9rRtGZ5Hh63SCObRaJEMAthJtARwMzqACeccwf892OBNkAO0AVoCkwJYjY5TWpyAmP7tOamptX482cr+e8Pv9MVSCIRLjZYb+ScW2ZmU81sJuADBppZDyDeOTfazE4BecAR4CXn3NZgZZMzS4zz8eK9zalZoQwvzVjH9v3HeLFrc8rE+7yOJiIBYOF6KiArK8vl5uZ6HSNqvDF3A09/9B1NalZgdK8s0pITvI4kIj+BmeU557LOtE0D0+S89Gybycj7WvL9jgPcMXwO6/IPeR1JREqZCkHOW8dGVXm7bxsOHy/izhFzyNtY4HUkESlFKgS5IM0zKjLxoWzKl4mj28vf8NlSXQwmEilUCHLBaqcmMfGhdjSqXo6Hxi1kzOz1XkcSkVKgQpCfpFJSPOP6tuH6BlV4+uPl/H3ySo1VEAlzKgT5yRLjfIzo3oKuV9Zi2NS1/GrCEopOnvI6loj8REEbhyCRKdYXwzO3NyE9OYEhX6+h4HAhL3ZtobEKImFIRwhy0cyMJzpexu9va8RXK3fRY/Q37DtS6HUsEblAKgQpNT3bZjK0awuWbNlPl5Fz2b5fU2iLhBMVgpSqm5pWY8wDrdi+/xh3Dp/Dml0HvY4kIudJhSClLvuSNN7p14bCk467Rs4lb6PWVRAJByoECYjGNcozcWDxALbur8zj65U7vY4kIj9ChSABk5FalvEDsqlXOZm+b+QxPm+L15FE5BxUCBJQ6SkJvNOvLW3qVuKX73/LSC22IxKyVAgScMkJsbzauxU3N63GXz5byZ8mreCUFtsRCTkamCZBkRDrY8i9zUlNiueVWespOFzIX+9qSpxPv5OIhAoVggRNTIzx9K2NSEtO4Nkpq9h7pJBh3VtQNl5/DUVCgX49k6AyMx65rj7P3N6E6avy6f6KRjWLhAoVgniiW+sMhndvwXdbD2hUs0iIUCGIZzo11qhmkVCiQhBPnT6qedEmjWoW8YoKQTzXuEZ5JgxsS7nE4mU5p6/K9zqSSFRSIUhIqJ2axPiBbamTlkSfMQv4cPFWryOJRB0VgoSMyimJvNO/DVmZFXnsncW8OktrNYsEkwpBQkq5xDjG3H8lNzSqwu8/Wc7fPtdazSLBokKQkJMY52N495Z0vbIWw6etZfCEpVqrWSQINERUQpIvxnjm9iakJSfw4tdr2HukkCFdm5MYp7WaRQJFRwgSssyMX3S8jKdvaciUFTvp+ep89h894XUskYgV1EIws8FmNtt/a3PattvNLNfM8szsF8HMJaGtd7s6/PPe5izatJd7XprLrgPHvI4kEpGCVghm1hDoDLQHugFDT3vKC8ANQFtgkJlVDlY2CX23NqvOq71bsangCHeOnMP63Ye9jiQScYJ5hNABmOyKbQRizaxcie1FQFkgEdgO/Ns8BmbWz38UkZufr8FL0aZD/XTe7tuGQ8eK6DJyDsu27vc6kkhECWYhpAL7Stw/5H/sB88Ci4HvgK2Anf4CzrlRzrks51xWenp64JJKyGpWqwLjB2aTEOvj3lHzmLNmt9eRRCJGMAthL5BS4n4FYA+AmWUADwOXAJnACaB3ELNJGLkkPZkJA7OpXiGR3q8t4NOl272OJBIRglkIM4GOAGZWBzjhnDvg35YIHAMOOedOAjuAA2d8FRGgavlE3uvfliY1yzNo3ELenLfR60giYS9oheCcWwZMNbOZwFvAQDPrYWZ9nHOrgDeB2WY2AygDvB2sbBKeKpSN580+rbn2sso89cEyXvhylUY1i1wEC9d/QFlZWS43N9frGBICTpw8xeAJS5mwcAs92tTm6Vsb4Yv5t6+gRAQwszznXNaZtmmksoS9OF8M/+jSlLTkeF6asY6CI4U8d3czEmI1qlnkQqgQJCKYGb/u3IDU5Hie+XQl+44U8lKPLJIT9Fdc5Hxp6gqJKP1yLuHZLs2Yt66ArqPmsfvQca8jiYQNFYJEnDtb1uTlni1ZvesgXUbOZXPBEa8jiYQFFYJEpGsvr8JbD7am4HAhd46Yw8oduopZ5MeoECRitaxdifcHtCXGjLtHzmXBhgKvI4mENBWCRLRLq6QwfmBb0lISuO+Vb5iyfKfXkURClgpBIl7NimUZPyCby6umMODNPN7L3ex1JJGQpEKQqFApKZ5xfduQfUkqT45fwohpazWqWeQ0KgSJGkkJsYzu1Ypbm1Xnr5+v5I+TVnDqlEpB5AcatSNRJT42hhfuuYJKSfGMnrWeHfuP8ezdzbRWswgqBIlCMTHG/9zSkBoVyvCnT1ew88AxXu6ZRcWkeK+jiXhKp4wkKpkZfXPqMqxbC5Zs3c8dI+awcY+W5ZTodsGFYGZJZqbja4kINzWtxrgHW7P3SCF3DJ/Dok17vY4k4pkfLQQzizGzbmY2ycx2ASuB7Wa23Mz+bmb1Ah9TJHCyMisxcWA2SQmxdH15Hl98t8PrSCKeOJ8jhKkUL235a6Cqc66Wc64y0B6YB/zVzO4LYEaRgKubnszEh7K5vGo5+r+Zx5jZ672OJBJ0P7pAjpnFOedOXOxzSpsWyJFAOFp4ksfeWcQXy3fyYPs6/FfnBsRosR2JIOdaIOdHjxB++B+9mf3TzM74LyPYZSASKGXifYy4ryW9szN5ZdZ6Bo1byLETJ72OJRIUF/Kl8kHgIzNLAjCzG8xsdmBiiXjHF2M8fWsjfntzQz7/bgfdXp5HweFCr2OJBNx5F4Jz7imKF76f5i+CJ4DBgQom4rU+7eswvFsLvtt2gDuGz2Zt/iGvI4kE1HkXgpldB/QFDgNpwKPOuZmBCiYSCm5sUo1xfdtw8FgRdwyfw5y1u72OJBIwF3LK6DfAb51zVwN3Ae+a2bUBSSUSQlrWrsgHg9pROSWBnqPn8878TV5HEgmICzlldK1zbpb/56XAjcAfAxVMJJTUqlSWCQ9lk10vjcETl/LMpys4qYnxJMKcz8C0s11ZtB247lzPEYkk5RLjeLVXFj3b1mbUjHUMeDOPw8eLvI4lUmrOa2CamT1iZhklHzSzeKCtmb0O9ApIOpEQE+uL4fe3NebpWxry1YqddBk5l+37j3odS6RUnE8hdAJOAm+b2Q9TVqwHVgNdgRecc2MCmFEk5PRuV4fRvVuxqeAItw2dzZIt+7yOJHLRfnSk8v95slkcxVcYHXXO7QtUqPOhkcoSCr7fcZAHxixgz+HjvHDPFXRqXM3rSCLndFEjlUu8yBPAAuAt4Hdm9oCZtTSzhFLKKRJ2LquawgeD2tGgWjkGvLmQYVPXaGlOCVsXctnpIxRfbtoH+BqoBTwFLDOzZefzAmY22Mxm+29tSjxew8ymlbgVmNljF5BNxDPpKQm83bcNtzSrzt8nf8/j7y7WdBcSli5kxbTvgbWu+Nef9cCHP2wws/I/9ofNrCHQmeJZUjOACUAWgHNuK3C1/3ktgeeB4ReQTcRTiXE+htx7BZdXTeEfX3zPuvzDjOrZkmrly3gdTeS8XcgRwi7gVTOrc/oG59z+8/jzHYDJrthGINbMyp3heS8Bj5xpwjwz62dmuWaWm5+ffwHRRQLPzBh0TT1e7pHF+t2HueXF2eRtLPA6lsh5u5BCWEbxEcW/zGyLmU0xs2cv4M+nAvtK3D/kf+x/mdnNwDrn3LdnegHn3CjnXJZzLis9Pf0C3lokeK5vWIV/PZRNUoKPrqO+4b0Fm72OJHJezvuUkXPubz/8bGaxwGVAkwt4r71AyVNLFYA9pz2nB/D6BbymSEiqXyWFDwe145G3F/HkhCUs336Ap25qQKxPy5hL6PpJfzudc0XOue+cc+9cwB+bCXQE8J92OuGcO3Dac64BpvyUTCKhpkLZeF7r3Yo+7eswZs4Ger02n72aRltCWNB+XXHOLaN41PNMii9dHWhmPcysD4CZVQIKnHP6FyMRI9YXw29vbsg/ujRjwfq93DZsNqt2HvQ6lsgZXdDAtFCigWkSbhZu2kv/sXkcOV7Es3dfQafGVb2OJFGoVAamicjFaZFRkY8fbk+9KikMeDOPv36+UjOmSkhRIYgEUdXyibzXvw3dWmcwYtpaer06X8tzSshQIYgEWUKsj2dub8Lf7mrK/A0F3PLiLE2OJyFBhSDikbuzajFhQDYAd42cy7sLtBKbeEuFIOKhJjXL8/Ej7WldpxK/mrCUX09cwvEizYMk3lAhiHisUlI8Y+6/kkHXXMLb8zdz98i5bN2nRXck+FQIIiHAF2P85w2X81KPlqzNP8wtL85i9prdXseSKKNCEAkhNzSqykcPtyM1KZ4eo79hyFerdWmqBI0KQSTE1E1P5oNB7fiPK2rw3JRV9H5tPrsPHfc6lkQBFYJICEpKiOXZu5vxlzuaMH99AZ3/OZNv1p0+F6RI6VIhiIQoM+PeKzP4YFA7khNi6fryPIZNXcMpnUKSAFEhiIS4BtXK8dEj7bmpafESnfePWaDRzRIQKgSRMJCcEMuQe6/gj//RmLnr9nDTkJnkbtBqbFK6VAgiYcLMuK9NbSYOzCY+NoZ7Rs3jpelrdQpJSo0KQSTMNK5RPLr5hkZV+PNnK+n12nx2HTzmdSyJACoEkTBULjGOYd1a8KfbG7NgQ/FVSFO/3+V1LAlzKgSRMGVmdG9dm48fbk9acgL3v7aAP3yyXHMhyU+mQhAJc/WrpPDBoHb0alub0bPWc8fwOazNP+R1LAlDKgSRCJAY5+N3tzXmlZ5ZbNt3lJuHzOLdBZsI1yVyxRsqBJEIcn3DKnz+eA7NMyrwqwlLefjtRew/esLrWBImVAgiEaZKuUTG9mnNk50uY/KyHZr2Qs6bCkEkAvlijIeursf7A9oS5zPufXkez3y6Ql84yzmpEEQiWPOMikx6tANdr8xg1Ix13DZ0Niu2H/A6loQoFYJIhEtKiOWZ25vwau8sdh8q5Nahsxg5fa3WWZB/o0IQiRLXXl6FL36ew3WXV+Evn62k66h5bC444nUsCSEqBJEoUikpnhH3teDZLs1Yvv0AnV6YwXsLNuvyVAFUCCJRx8y4s2VNPn+8A01qlufJCUvoNzaP/INalS3aqRBEolTNimUZ92AbnrqpAdO/z6fj89P5cPFWHS1EsaAWgpkNNrPZ/lub07ZdYWbTzOwbMxtvZonBzCYSjWJijAc71GXSo+3JSE3isXcW039snmZPjVJBKwQzawh0BtoD3YChpz1lFNDLOdca+BrIDFY2kWhXv0oKEwa05dc3Xs60Vfl0fH4GHyzS0UK0CeYRQgdgsiu2EYg1s3IAZpYJHAF+ZWYzgBTn3MrTX8DM+plZrpnl5ufnBzG6SOSL9cXQ/6pL+PTRDtRJS+LxdxfT9408dh3Q0UK0CGYhpAL7Stw/5H8MoBrQBhgJXAtcY2bXn/4CzrlRzrks51xWenp6gOOKRKd6lZMZPyCb33RuwMzV+fzs+RlMXLhFRwtRIJiFsBdIKXG/AvDDBCvHgLXOuSXOuSJgEtA8iNlEpARfjNE3py6fPtaBepWTeeK9b3nw9Vx27NfRQiQLZiHMBDoCmFkd4IRz7ocx9CuAVDOr679/FbAkiNlE5AwuSU/mvf5teeqmBsxas5vrn5vO2LkbtI5zhApaITjnlgFTzWwm8BYw0Mx6mFkf59wxoAcwzszmABucc5ODlU1Ezs7nvxJp8uM5NKtVnt9++B1dXprLqp0HvY4mpczC9bxgVlaWy83N9TqGSFRxzjFx4Vb+MGk5h48XMfCqS3jomnokxvm8jibnyczynHNZZ9qmgWkict5+GOX81RNXcXPT6gz5eg2dh2i9hUihQhCRC5aanMDz91zB6w9cSWHRKe4ZNY/BE5aw/4hWZwtnKgQR+cmuujSdL36eQ7+curyXu5nrnpvOx99u0yWqYUqFICIXpWx8LP/VuQEfPdyeauUTeeTtRfQYPZ+1+Ye8jiYXSIUgIqWicY3yfDCoHb+/rRHfbtlHpxdm8PfJKzlaqGU7w4UKQURKjS/G6Nk2k69/cTW3NKvOsKlruf656UxZvtPraHIeVAgiUurSUxJ47u4reLdfG5ISfPR9I5c+YxZohbYQp0IQkYBpXTeVSY924DedGzBv3R6uf246Q75azfEinUYKRSoEEQmoOF8MfXPq8uUvruL6BlV4bsoqbnh+Bl+t2KmrkUKMCkFEgqJa+TIM696CsX2uxBdj9Hk9l96vLWDNLk2BESpUCCISVB3qp/P54zn8980NWbhpL51emMnvP17O/qMa1OY1FYKIBF2cL4YH2tdh2i+v5u5WtXhtznqu+cc03vpmIyc1k6pnVAgi4pnU5ASeub0JnzzSnnqVk/nNv5Zx84uzmKe5kTyhQhARzzWqXp53+7VhWLcWHDh6gntHzWPQWwt1mWqQxXodQEQEimdSvalpNa5rUJlRM9YxfNoapizfSa/s2jx8TX3Kl43zOmLE0xGCiISUxDgfj15Xn6m/vJrbrqjOK7PWk/P3qYyetZ7ColNex4toKgQRCUnVypfh712aMemRDjStWZ4/fLKc65+bzqQl2zV+IUBUCCIS0hpWL8fYPq15/YErKRvvY9C4hdwxYg65Gwq8jhZxVAgiEhauujSdSY924G93NmXbvqPcNXIuA8bmsX73Ya+jRQx9qSwiYcMXY9zdqhY3N6vGKzPX89L0tXy5Yif3tKrFo9fVp0q5RK8jhjUL13NxWVlZLjc31+sYIuKh/IPHGfLVat6ev4lYn9ErO5MBOZdQMSne62ghy8zynHNZZ9ymQhCRcLdpzxGe/3IVHyzeSnJ8LP1y6vJA+zokJegkyOlUCCISFVbuOMCzX6xiyvKdpCbFM+iaenRvk0FCrM/raCFDhSAiUWXhpr38/fPvmbtuDzUqlOGx6+tzR/MaxPp0Hc25CkH/dUQk4rTIqMi4vq15s09rUpPjeXL8Ejq+MIMPF2/V5HnnoEIQkYhkZrSvn8aHg9ox8r4WxMXE8Ng7i+n4/HQVw1moEEQkopkZnRpX47PHOjCsWwt8McZj7yzmhhdm8NG321QMJagQRCQqxMQUT573+WM5DO3WHAMefXsRnV6YwcffbuOUiiG4hWBmg81stv/W5rRtj5jZcjOb5r9dFsxsIhIdYmKMm5tWZ/LjxcUA8Mjbi7jhhRl8siS6iyFoF+maWUOgM9AeyAAmACW/6W4J9HbOzQ9WJhGJXj8Uw42Nq/Hp0u3886vVPDxuEfUrr2bQNfW4uWm1qLsqKWiXnZpZfyDNOfcn//3FQI5z7oD//lJgNVAZmOSc+/MZXqMf0A8gIyOj5caNG4OSXUQi38lTjklLtzP069Ws2nmIWpXKMOCqS7izRU0S4yJnHEOoXHaaCuwrcf+Q/7EfvA/0B64F2pvZrae/gHNulHMuyzmXlZ6eHsisIhJlfDHGrc2q8/ljOYzq0ZJKZeP5zb+WkfO3qbwycx2Hjxd5HTHgglkIe4GUEvcrAHsAzMyA551z+c65QmAS0DSI2UREgOJTSR0bVeWDQe1468HW1KuczB8nraDdX7/mn1+uZv+RE15HDJhgFsJMoCOAmdUBTvxwughIBlaaWTl/OVwP6LsEEfGMmdGuXhrj+rZh4kPZZNWuyPNfriL7L1/x589WsOvgMa8jlrqgTl1hZr+luBR8wBNAfSDeOTfazLoDjwOFwBTn3NPnei1NXSEiwbZi+wFGTFvLJ0u2EeuL4a6WNXmwfR3qpid7He28aS4jEZFStGH3YV6asZYJC7dy4uQprm9QhX45dcmqXZHikxyhS4UgIhIA+QePM3buBt6Yt5F9R07QPKMC/TrUpWOjqvhiQrMYVAgiIgF0pLCI8XlbeGXmejYVHCGjUlke7FCHu1rWpGx8aK3JoEIQEQmCk6ccX3y3g5dmrGPx5n1UKBtHjza16dk2k/SUBK/jASoEEZGgcs6Ru3Evo2as48sVO4mLieGWZtW5v10mjWuU9zTbuQohtI5lREQigJnRKrMSrTIrsTb/EGNmb2DCwi1MWLiFVpkVub9dHTo2rBJyU2PoCEFEJAj2Hz3B+7mbeX3uBjYXHKV6+UR6tM3k3la1qJgUH7QcOmUkIhIiTp5yfLViJ2PmbGDO2j0kxMZwe/Ma9G6XyeVVywX8/VUIIiIhaOWOA7w+ZwMTF27leNEp2tZNpXe7TK67vHLATiepEEREQtjew4W8s2AzY+duYNv+Y1Qrn0i3KzO458paVE5JLNX3UiGIiISBopOn+HLFLt76ZiMzV+8mNsa4oXFV7mtdmzZ1K5XKKGhdZSQiEgZifTF0alyVTo2rsn73Yd6at5H387Ywacl26lVO5r7WGdzRsiblEuMC8v46QhARCWHHTpzk42+38eY3m/h28z7KxPl44meX0jen7k96PR0hiIiEqcQ4H12yatElqxZLt+znzXkbqVGxTEDeS4UgIhImmtQsz1/vCtzaYaE1TE5ERDyjQhAREUCFICIifioEEREBVAgiIuKnQhAREUCFICIifioEEREBwnjqCjPLBzb+xD+eBuwuxTjhQPscHbTP0eFi9rm2cy79TBvCthAuhpnlnm0uj0ilfY4O2ufoEKh91ikjEREBVAgiIuIXrYUwyusAHtA+Rwftc3QIyD5H5XcIIiLy76L1CEFERE6jQhARESAKCsHMBpvZbP+tzWnbmpvZTP/td15lLG0/ss+PmNlyM5vmv13mVc7SZGY5ZjbjDI9H5GcM59zniPuMzSzWzMb4P8f5ZnbLads7mtlc/9/5/l7lLE3nsc/PmtnCEp9z+Yt+U+dcxN6AhsAMwIDaQO5p2+cBl/p/ngw09zpzEPZ5DHCl1zlLeZ9/BSwB5p1hW8R9xuexz5H4GfcChvp/TgM2ltgWC6wAUoE4YBFQ2evMgdxn/2PTSns/I/0IoQMw2RXbCMSaWTkAM0sAKjnnVvmf+5n/+eHurPvs1xIYbGazzOzX3kQsdWuAO09/MII/YzjLPvtF4mc8AfiN/+dTp227BNjinNvjnDsBTAfaEP7Ous9mZsBlwEv+z/n+0njDSC+EVGBfifuH/I/9sG3/WbaFs3PtM8D7QH/gWqC9md0avGiB4ZybAJw4w6ZI/YzPtc8QmZ/xIefcfjNLAcYDvy2x+cf+zoelH9nnssBwoDvQCRhkZs0u9j0jvRD2Aikl7lcA9pzHtnB21v3y/1bxvHMu3zlXCEwCArdit/ci9TM+q0j+jM2sBvAlMM4590aJTRH7OZ9jn48BzznnjjjnDgFfA00u9v0ivRBmAh0BzKwOcMI5dwDAOXcU2G9mdf3/iG4EZnmWtPScdZ+BZGClmZXz7/P1wHxvYgZeBH/G5xKRn7GZVQO+AH7jnHvltM2rgdpmVsHM4oEc4JtgZyxtP7LP9YBv/F88x1F8KjTvYt8z9mJfIJQ555aZ2VQzmwn4gIFm1gOId86NBh4G3qD4C9gvnXMLPYxbKn5sn83sSeAroBCY4pz7wsu8gRDpn/GZRMFnPJji00BPmdlT/se+AhY65yaZ2S+ATyn+JXeYc26nRzlL04/t82vAXKAIeM05t+Ji31AjlUVEBIj8U0YiInKeVAgiIgKoEERExE+FICIigApBRET8VAgiIgKoEERExE+FIFJK/AMCf+b/+Y9m9qLXmUQuRESPVBYJsv8Bfm9mlYHmQNhPKifRRSOVRUqRmU2neD6hq51zB73OI3IhdMpIpJSYWROgGlCoMpBwpEIQKQX+mSnfAm4DDplZJ48jiVwwFYLIRTKzssBE4Bf+GSf/QPH3CSJhRd8hiIgIoCMEERHxUyGIiAigQhARET8VgoiIACoEERHxUyGIiAigQhAREb//Bz0NgCQ/pWyoAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "delta = 0.0001\n",
    "x = np.linspace(delta, 2.5, 1000)\n",
    "mH = np.sqrt(2) * x\n",
    "I_0 = iv(0, 2*mH)\n",
    "I_1 = iv(1, 2*mH)\n",
    "I_2 = iv(2, 2*mH)\n",
    "eta = 2 / mH * I_2 / I_1\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(x, eta)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$\\eta(x)$')\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.13"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}